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L  Introduction 


Over  the  past  40  years,  studies  have  been  made  [1-4]  of  the  process  of  ob¬ 
taining  images  of  the  reflectivity  or  density  of  target  areas  that  are  rotating 
and  translating  with  respect  to  a  sensor  such  as  a  monostatic  or  bistatic  ra¬ 
dar,  a  sonar,  or  an  x-ray  CAT  scanner.  Ausherman  et  al  [5]  have  written  an 
excellent  review  of  the  work  done  in  this  area.  Target  areas  that  rotate  and 


translate  relative  to  a  radar  include,  for  example,  planet  surfaces  observed 
from  a  satellite,  ground  terrain  observed  from  airborne  platforms,  subsur¬ 
face  objects  and  voids  observed  from  moving  vehicles,  and  people  scanned 
by  a  rotating  x-ray  system.  Although  the  processing  described  is  applicable 


to  other  systems,  this  report  treats  the  topic  from  the  point  of  view  of  a  syn¬ 
thetic  aperture  radar  (SAR), 


In  the  SAR  application,  as  the  aspect  angle  between  the  sensor  and  a  target 
changes  with  time,  the  sensor  collects  a  sequence  of  signal  records.  After 
data  are  collected  for  Ts  seconds,  the  aspect  angle  has  changed  by  &s  de¬ 
grees.  These  received  signals  are  then  coherently  processed  to  produce  the 
reflectivity  profile  of  the  target  area.  Down-range  resolution  into  range  bins 
is  determined  primarily  by  the  bandwidth  of  the  sensor.  Cross-range  resolu¬ 
tion  is  obtained  primarily  by  coherent  processing  of  the  received  signals  so 
that  a  very  wide  aperture  is  simulated:  an  aperture  that  is  0S  degrees  wide. 


Implementation  and  study  of  image  formation  processing  have  been  limited 
in  two  wavs.  First,  the  image  formation  processing  has  assumed  that  targets 
are  isotropic  point  scatterers.  although  many  targets  are  anisotropic  reso¬ 
nant  seatterers.  Treatment  of  resonant  scattering  becomes  important  when 
the  sensor  spectrum  covers  the  Rayleigh,  resonant,  and  optical  regions  of  a 
family  of  targets.  For  example,  discrimination  between  scatterers  can  be 
based  on  the  unique  signature  of  each  target.  But  in  order  for  the  discrimi¬ 
nation  to  work  in  the  context  of  microwave  reflectivity  imaging,  the  infor¬ 
mation  in  the  signature  must  be  preserved  during  the  image  formation 
process. 


The  second  limitation  in  previous  work  is  that  systems  with  relatively  nar¬ 
row  bandwidth  have  been  assumed.  For  example,  the  compressed  pulse 
width  of  the  sensor  is  assumed  to  be  at  least  several  and  usually  many 
cycles  of  an  rf  carrier  frequency.  So  a  single  range  bin  is  derived  from  sev¬ 
eral  cycles  of  the  carrier.  Newer  UWB  (ultra-wide-bandwidth)  sensors, 
however,  have  made  it  possible  to  make  the  image  range-bin  size  roughly 
half  the  wavelength  of  the  highest  frequency  in  the  sensor’s  spectrum.  The 
bin  size  is  much  smaller  (1/10  or  less)  than  the  wavelength  of  the  lowest 
frequency  in  the  spectrum.  This  wide  bandwidth  exacerbates  range  walk 
and  wavefront  curvature  errors  to  the  point  where  conventional  fast  Fourier 
transform  (FFT)  based  processing  must  be  restricted  to  very  small  patches 
within  an  image  area. 
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These  two  limitations  are  not  independent.  It  is  bandwidth  that  allows  target 
discrimination  based  on  signature  analysis,  and  it  is  bandwidth  that  makes 
image  formation  more  difficult.  The  purpose  of  this  report  is  to  examine 
image  formation  processing  that  preserves  resonant  target  signatures  and  to 
present  an  efficient  method  of  solving  the  microwave  reflectivity  imaging 
problem  for  UWB  signals  and  resonant  targets. 


2,  Background 


2.1  Terminology 

SAR  systems  depend  upon  collecting  data  coherently  along  a  path.  This 
path  is  referred  to  as  the  “synthetic  aperture.”  Figure  1  shows  a  sketch  of 
the  scenario  projected  onto  a  plane.  An  aircraft  carrying  the  SAR  system 
flies  at  some  elevation  h  above  the  ground,  over  a  distance  JL_.  to  form  the 

b' 

synthetic  aperture.  The  image  area  grid  is  referenced  to  the  center  of  the  ap¬ 
erture.  Range  is  marked  off  by  the  parameter  i.  The  parameter  k  marks  off 
azimuthal  (bearing)  lines.  The  sample  points  in  the  aperture  are  marked  off 
by  the  parameter  j.  The  distance  between  the  ;'th  point  in  the  aperture  and  the 
(i,k)th  position  in  the  image  area  is  denoted  by  d- Both  the  azimuthal 
lines  and  the  range  lines  are  referenced  to  the  center  of  the  aperture.*  Table 
1  summarizes  the  nomenclature  used  in  this  report. 

(Although  nearly  all  operational  SAR  systems  use  a  rectangular  grid,  the 
advantages  of  a  polar  grid  are  that  for  a  ringing  target,  the  grid  provides 
equally  spaced  samples  of  the  ring-down,  and  these  samples  are  optimally 
focused  and  on  a  single  bearing  line.  These  features  make  it  easy  for 
postprocessors  to  look  for  specific  ringing  signatures.) 


2.2  Approaches  to  Focusing 

A  “Doppler”  paradigm  is  often  used  in  discussions  of  SAR  processing,  or 
SAR  focusing.  Referring  to  figure  1,  assume  that  the  radar  platform  is  mov¬ 
ing  at  velocity  v  as  it  collects  data  along  the  aperture.  Therefore,  the  data 
collection  points  are  spaced  equally  in  time,  occurring  at  a  rate  governed  by 
v  and  the  PRF.  (This  sampling  across  the  aperture  is  sometimes  referred  to 


Figure  1.  Basic  flight 
path  and  image  area 
geometry. 


Flight  path 


*  This  grid  is  convenient  for  postprocessing.  An  x-y  grid  can  also  be  computed  with  the  techniques  described. 
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Table  1.  Analysis  nomenclature. 


Symbol 

~ 


Caret— e.g.,  s(t) 
Subscript  i,k,j 


c 

a  ~  h±\,k,j~  hkj 

PRF 

u 

UWB 


Meaning 

Length  (meters)  of  the  synthetic  aperture. 

Denotes  a  received  signal  amplitude  (volts)  as  a  function  of  time  (seconds  from  the  leading 
edge  of  the  transmitted  pulse),  at  the  position  in  the  aperture.  The  analog-to-digital  outputs 
represent  this  signal. 

Denotes  that  the  function  is  an  estimate.  The  example  denotes  an  estimate  for  a  received 
signal. 

Particularizes  the  function  to  be  at  particular  geometries,  like  the  j,h  position  in  the  synthetic 
aperture,  or  the  £th  bearing  line,  or  the  ( i,k)th  position  in  the  image  area. 

Denotes  the  distance  (meters)  between  the /h  position  of  the  radar  and  the  (i,k)lh  position  in  the 
area  to  be  imaged. 

Denotes  the  round-trip  time  (seconds)  for  the  radio  energy  to  travel  from  the/*'  position  in  the 
synthetic  aperture,  to  the  (i,k)ih  position  in  the  image  area,  and  back. 

The  speed  of  light. 

Time  shift  between  range  bin  i  and  range  bin  i  +  1  at  the  center  of  the  aperture  (where  j  =  0). 
Pulse  repetition  frequency  of  the  radar  in  hertz. 

Relative  bandwidth  u  -  (Fh:  —  F[0)IFq  and  Fq  ~  ( F y.  +  F-lo)j1. 

Ultra-wide  bandwidth,  where  1. 


fa  (') 


Wavelength  in  meters. 

Impulse  response  of  target  at  (i,k). 


as  “slow  time.”)  If  the  radar  is  operating  at  a  frequency  /q,  then  the  echo 
from  a  target  in  the  image  area  will  have  a  Doppler  shift  profile  as  the  plat¬ 
form  moves  through  the  aperture.  The  shift  will  be  upward  while  the  plat¬ 
form  approaches,  it  will  drop  to  zero  as  the  platform  moves  to  a  position 
broadside  to  the  target,  and  it  will  be  downward  as  the  platform  recedes. 
Every  target  position  will  have  a  unique  Doppler  profile.  Since  every  posi¬ 
tion  in  the  image  has  a  distinct  Doppler  profile,  one  can  form  an  image  by 
simply  assigning  to  each  pixel  a  filter  matched  to  the  Doppler  profile  ex¬ 
pected  for  a  target  at  the  location  represented  by  that  pixel. 

The  classic  “polar  formatting”  [6]  approach  to  computing  a  SAR  image  ad¬ 
justs  (time  shifts)  the  data  at  each  aperture  point  so  that  the  new  data  set  ap¬ 
pears  as  if  the  platform  had  moved  in  a  short  circular  path,  with  the  circle  at 
the  center  of  the  image  area.  Once  this  formatting  is  done,  a  target  at  the 
center  point  has  the  unique  Doppler  profile  of  zero  over  the  entire  aperture. 
Other  points  have  other  unique  Doppler  profiles. 

A  Fourier  transform  is  typically  used  to  recover  the  image  from  the  Doppler 
profiles.  This  approach  works  well  as  long  as  the  circular  path  is  suf¬ 
ficiently  short,  the  image  area  sufficiently  small,  the  signal  bandwidth 
sufficiently  small,  and  the  distance  from  the  radar  to  the  image  center  suffi¬ 
ciently  long.  This  report  addresses  the  case  where  none  of  these  restrictions 
apply. 
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Although  the  Doppler  paradigm  has  helped  countless  people  to  visualize 
SAR  focusing,  it  is  not  necessary,  nor  always  helpful,  in  solving  the  SAR 
focusing  problem.  For  example,  Doppler  shift,  which  is  defined  as  2v/A, 


works  fine  when  A  varies  by  a 


few'  percentage  points,  but  it  becomes  a 


stumbling  block  in  the  UWB  case  where 
although  Doppler  has  proven  to  be  a  ver 


A  can  vary  by  10  or  100  to  1.  Thus, 
y  convenient  narrow-band  concept. 


its  convenience  breaks  down  at  wide  relative  bandwidths. 


The  focusing  problem  can  also,  however,  be  looked  at  as  a  stationary  array 
of  N  antennas  whose  outputs  are  digitally  stored  and  combined  in  a  com¬ 
puter  to  form  beams.  We  believe  that  the  focusing  problem  is  easier  to  visu¬ 
alize  and  solve  within  this  stationary  array  paradigm  when  neither  the  ge¬ 
ometry  nor  the  bandwidth  are  restricted.  In  this  paradigm  one  can  see  that 
equation  (1)  coherently  focuses  the  SAR  data  by  summing  across  the  array: 


f  i,k  (0  =  2  s}  (TiikJ  + 1)  for  t  *  0  .  (1) 

j 

Both  the  Doppler  and  stationary  array  approaches  can  be  seen  to  be  identi¬ 
cal  at  the  center  point  of  the  “polar  format”  patch.  Notice  that  regardless  of 
carrier  frequency,  the  center  point  of  the  polar  formatted  data  is  always  ana¬ 
lyzed  by  the  dc  term  of  the  Fourier  transform,  which  is  simply  a  summation 
of  the  data  points.  The  summation  shown  in  equation  (3.)  is  identical  to  find¬ 
ing  the  dc  term  of  the  polar  formatted  SAR  data.  But  instead  of  formatting 
once  and  then  finding  many  “Doppler”  profiles  via  an  FFT,  equation  (!) 
formats  and  “sums”  the  data  many  times:  each  formatting  makes  a  different 
pixel  the  center,  and  then  the  dc  term  summation  is  used  to  calculate  the 
value  for  that  pixel.  The  result  is  that  optimally  focused  beams  are 
formed— beams  taking  full  advantage  of  both  the  entire  aperture  and  the 
entire  signal  bandwidth.  The  focusing  is  truly  frequency  independent. 

How  does  the  beam  width  compare  with  that  obtained  from  “conventional” 
(far-field  narrow-band)  antenna  theory?  The  rule-of-thumb  half-power 
beamwidth  of  a  line  array  is  approximately  XfL,  where  L  is  the  length  of  the 
array.  The  beams  formed  by  the  processing  described  above  follow  this  rule 
and  have  a  width  proportional  to  A.  Low  frequencies  have  wide  beams,  and 
high  frequencies  have  narrow'  beams.  This  fact  has  interesting  conse¬ 
quences  in  the  time  domain  as  a  source  moves  through  a  beam.  Although 
the  impulse  response  at  the  center  of  the  beam  is  a  narrow  pulse,  as  one 
moves  away  from  the  center,  the  impulse  response  broadens.  This  time-do¬ 
main  broadening  occurs  because  more  and  more  high-frequency  energy  is 
lost  as  the  source  moves  out  of  the  narrowing  high-frequency  beam. 

The  switch  from  the  Doppler  to  the  stationary  array  approach  may  lead  to 
thinking  in  “phased  array”  terms.  Such  an  approach  is  a  mistake  when 
bandwidth  and/or  geometry  are  not  restricted.  Whenever  the  geometry  is 
not  restricted  to  the  far  field  of  an  aperture,  plane-wave  simplifications 


9 


break  down.  This  breakdown  invalidates  simple  phase  steering.  Since  a 
Fourier  transform  forms  beams  by  simple  phase  steering,  Fourier  tech¬ 
niques  becomes  less  and  less  useful  as  targets  move  into  the  near  field.  In  a 
UWB  system,  phase  becomes  meaningless  with  regard  to  defining  element 
positions  or  the  beam-forming  network.  To  speak  of  shifting  one  element 
180°  with  respect  to  another  element  implies  a  fixed  A,  If,  for  example,  A 
changed  by  2  to  1,  then  a  delay  line  that  provided  180°  at  one  frequency 
would  provide  360°  at  the  other.  Yet  delay  lines  are  precisely  the  element 
needed  to  build  a  wide-bandwidth  “phased-array”  antenna.  When  A  varies 
several  octaves,  the  best  parameters  to  use  in  the  equations  defining  the  an¬ 
tenna  are  the  time  and  distance:  the  time  shift  of  the  delay  lines  that  form 
the  beam-forming  network,  and  distance  between  antenna  elements.  Thus  it 
is  best  to  think  in  “timed-array”  terms.  This  time-based  framework  results 
in  derivations  that  are  frequency  and  geometry  independent. 


Null  Steering,  Pattern  Forming,  and  Grating  Lobes 

Generally,  an  antenna  designer  would  say  that  the  most  critical  aspect  of 
making  desirable  antenna  patterns  is  forming  nulls.  A  simple  classic  case  is 
spacing  two  elements  at  90°  and  phasing  them  by  90”  to  form  an  endfire 
beam  in  one  direction  and  a  null  in  the  opposite  direction.  In  the  simplest 
case  (where  the  element  weighting  is  +1  and  —1),  to  form  a  UWB  null  one 
could  time-steer  the  array  to  where  the  null  should  be,  and  then  invert  half 
the  elements  before  summing.  Of  course,  using  other  weighting  factors  al¬ 
lows  more  freedom.  An  impulse  signal  6(t)  coming  from  a  direction  other 
than  the  null  would  produce  in  the  receiver  some  array-induced  waveform. 
If  we  ignore  bandwidth  effects  from  each  element,  that  waveform  would  be 
a  function  of  the  +  and  -  ’weighting  and  element  spacing.  By  definition,  that 
waveform  is  the  antenna-arrav  impulse  response  for  that  beam  angle. 
Matched  filtering  to  that  waveform  forms  a  main  lobe  in  that  direction. 

Grating  lobes  are  customarily  defined  as  lobes  whose  gain  is  equal  to  the 
main  beam.  However,  this  definition  leads  to  confusion  in  the  UWB  case. 
Are  there  grating  lobes?  That  depends.  If  the  response  to  the  UWB  signal  is 
seen  as  an  impulse,  then  the  answer  is  no.  Since  the  peak  response  on  the 
main  lobe  is  higher  than  the  response  at  any  other  angle,  the  definition  fails. 
On  the  other  hand,  the  antenna  is  a  linear  system.  It  behaves  just  like  an 
identical  array  operating  at  a  single  frequency.  So,  if  the  response  is  to  a  cw 
signal,  the  answer  can  be  yes.  If  the  elements  are  physically  spaced  more 
than  A/2  apart  at  the  highest  frequency  component  in  the  UWB  waveform, 
then  there  certainly  will  be  a  grating  lobe  at  that  frequency  component.  All 
the  insight  gained  from  cw  antenna  analysis  holds  and  remains  useful.  One 
must,  however,  be  careful  when  applying  terms  like  grating  lobes  that  pre¬ 
suppose  a  narrow-band  signal.  For  SAR,  the  element  spacing  (velocity/ 
PRF)  needed  will  be  a  function  of  what  sidelobes  are  permissible  at  the 
various  frequency  components  in  the  UWB  waveform. 


2.4  Target  Resonance  Effects 

Historically,  the  relative  bandwidth  of  radars  has  been  sufficiently  small 
that  a  target’s  echo  is  adequately  modeled  by  a  single  number,  a,  the  radar 
cross  section  (RCS),  usually  given  in  square  meters.  When  ji  ss  0.5,  how¬ 
ever,  a  single  number  may  no  longer  be  adequate — o  is  a  function  of  fre¬ 
quency.  For  example,  figure  2  is  a  plot  of  the  RCS  of  a  sphere.  In  addition 
to  the  magnitude  characteristic  plotted,  there  is  also  a  phase  characteristic. 
These  two  frequency-domain  characteristics  can  also  be  represented  in  the 
time  domain  by  a  ringing  or  resonant  response.  In  either  case — time  domain 
or  frequency  domain — the  waveform  in  the  plots  can  be  referred  to  as  the 
impulse  response  of  the  target. 


Since  historical  SAR  systems  have  not  been  designed  to  respond  to  target 
ringing,  no  attention  was  paid  to  preserving  the  ringing  information.  This 
report  describes  and  analyzes  a  procedure  to  focus  a  large  array,  over  ultra¬ 
wide  bandwidths,  in  such  a  way  as  to  preserve  the  resonant  response  of  tar¬ 
gets,  even  when  they  are  in  the  near  field. 


Figure  2.  Radar  cross 
section  of  a  sphere. 
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3.  Fundamental  Time-Based  Physical  Array  Approach 

Consider  an  aperture  looking  at  completely  empty  space  except  for  an  iso¬ 
tropic  scauerer  at  position  (i,k).  Also  assume  that  an  ideal  impulse  3(;)  is 
broadcast.  It  is  desirable  to  take  advantage  of  all  the  echo  energy  across  the 
aperture.  To  do  that,  one  must  sum  the  energy  from  all  collection  points 
along  the  synthetic  aperture.  However,  this  summation  must  be  done  so  that 
the  energy  adds  coherently  at  all  frequencies.  Frequency-independent  add¬ 
ing  is  accomplished  as 

fi* (0 *  £ sj(Tixj  + ?) wj for t~°  ■  a ) 


Her would  be  the  impulse  response  of  the  target  located  at  position 
(i,k).  Note  that  the  Ti>kj  term  time  shifts  the  received  signals  s{  so  that  the 
target  impulse  response  starts  at  t  =  0  at  all  points  in  the  aperture.  The  Wj 
term  is  simply  an  amplitude  tapper  across  the  aperture. 


The  above  discussion  assumes  that  a  target’s  impulse  response  is  suffi¬ 
ciently  similar  at  all  points  along  the  aperture  for  equation  (1)  to  be  consid¬ 
ered  true.  The  response  of  a  vertical  dipole,  for  example,  does  not  change 
depending  on  where  the  radar  is  positioned  in  the  aperture.  It  is  therefore  an. 
isotropic  target. 

Suppose  we  relax  the  isotropic  requirement,  and  suppose  that  a  complex 
scatterer  is  at  position  (i,k)  in  the  image  area.  A  horizontal  dipole,  for  ex¬ 
ample.  is  an  anisotropic  target.  We  could  take  advantage  of  the  anisotropic 
behavior  to  extend  the  detection  and  target  recognition  performance  of  the 
radar  by  rewriting  equation  (1)  as 

/,-,*( 0  -  Jy  (TiXi +  0 for  * *  /2\ 

where  y  ft)  =  sfi)  ©  xfj)  and  the  convolution  step  represents  filtering.  In 
this  case,  a  set  of  filters  X-(w)  would  be  needed,  each  one  matched  to  the 
target  response  at  the  bearing  of  the  ftt  position  in  the  aperture.  In  order  to 
simplify  the  rest  of  this  report,  we  use  equation  (1)  as  the  fundamental 
equation.  Nonetheless,  one  roust  recognize  that  real  targets  are  complex 
anisotropic  polarimetric  scatterers;  this  fact  should  not  be  ignored  for  large 
arrays. 
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4*  Short  Impulse  Response  Approximation 

A  typical  2-D  SAR  image  is  simply  the  echo  magnitude  mapped  to  inten¬ 
sity.  Equation  (1)  expands  the  typical  2-D  image  into  a  3-D  image  with  the 
target  ringing  along  the  third  dimension.  Given  the  heavy  computational 
load  of  typical  SAR  processing,  if  it  is  required  that  an  f(t)  (rather  than  a 
single  value)  he  computed  for  every  pixel,  then  a  massive  computer  would 
be  needed.  The  impact  of  this  computational  load  is  even  worse  when  one 
considers  that  the  mass  of  resulting  data  must  be  analyzed  in  the  target- 
detection/identification  phase.  This  section  defines  an  approximation  that 
reduces  the  problem  back  to  a  2-D  case  and  presents  an  error  analysis  of  the 
approximation. 

4.1  Theory 

The  problem  is  that  we  have  added  a  third  dimension  to  measure  target 
ringing.  Note,  however,  that  if  there  was  only  one  aperture  position,  then 
ringing  of  the  target  would  appear  in  pixels  only  behind  the  target.  Even 
with  an  aperture  of  many  positions,  ringing  will  appear  behind  the  target. 
But  since  the  geometry  is  not  constrained,  and  near-field  operation  is  pre¬ 
sumed,  the  antenna  beam  defocuses  behind  the  target. 

Equation  (1)  will  perfectly  focus  ringing  of  unbounded  duration.  For  practi¬ 
cal  purposes,  however,  the  target  rings  only  for  a  finite  duration,  say  M 
range  bins.  The  question  is,  given  that  the  ringing  is  finite  in  duration,  can 
the  ring  information  can  be  retained  in  a  2-D  image?  In  other  words,  sup¬ 
pose  that  instead  of  calculating  a  time  series  f;  k(t)  for  each  pixel  in  the  im¬ 
age,  only  one  value  is  calculated,  for  example,  fi  k(t),  where  r  is  fixed  for 
the  image.  If  the  ringing  can  be  retained  in  this  2-D  image,  then  computa¬ 
tional  load  is  greatly  reduced.  The  aim  of  this  section  is  to  describe  and 
quantify  the  error  bounds  when  a  2-D  image  /■  k{x)  is  used. 

First,  define  a  to  be  the  round-trip  time  it  takes  the  radar  pulse  to  traverse 
one  range  bin  on.  the  polar  grid;  that  is, 

«  “  (M,kj  ~  WVfc>  and  /  -  0  •  (3) 

Since  the  grid  for  the  image  has  been  defined  to  be  referenced  to  the  center 
of  the  aperture  (the  i  =  0  position),  the  i  parameter  can  be  thought  of  as  a 
quantized  time  t  parameter.  While  t  is  in  units  of  seconds,  i  is  in  units  of 
range  bins,  and  they  are  related  by  a — one  range  bin  equals  a  seconds. 
Now  we  can  write 


sj(TiXj  +  ma)  m 


for 


f  case  1:  j  «  0,  Vm,  Vi,  Vk  | 
|  case  2:  m  -  0,  V),  Vi,  ¥k  j 


(4) 


%  0; 
f  >1l 
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We  make  an  approximation  to  equation  (4),  generalizing  it  to  include  all 
aperture  points  over  a  limited  range  of  m  to  arrive  at 

Sj(Ti  k  j  +  ma)  »  Sj(TUmiktj)  for  \fi,Vk,Vj,m  =  .  (5) 

The  following  may  be  said  of  the  approximation  in  equation  (5). 

It  is  perfect  at  the  center  of  the  aperture  (j  =  0)  regardless  of  m. 

It  is  perfect  at  m  =  0  regardless  of  j. 

It  gets  worse  as  m  deviates  further  from  zero. 

It  is  worst  at  the  end  points  of  the  aperture. 

A  solution  is  desired  to  equation  (1)  in  the  form  offi>k(r),  where  r  is  a  con¬ 
stant.  In  a  practical  system,  rwill  be  mapped  to  discrete  range  bins,  so  let 

t  —  net  .  (6) 

Substituting  equations  (5)  and  (6)  into  (1),  we  define  a  2-D  “image” /(/,&)  as 


j 

-  fi-n,k(na)  “  ^Sj(T'-n,k,j  +  na)  K 


(7) 


If  one  thinks  of  a  ringing  target,  then  equation  (7)  can  be  described  as  al¬ 
lowing  the  point  of  perfect  focus  to  be  adjusted  to  any  depth  n  in  the  ring. 
For  example,  if  n  =  0,  then  the  perfect  focus  point  would  be  at  the  leading 
edge  (the  first  sample)  of  the  ringing  response.  If  n  *  0,  say  n  =  3,  then  the 
third  sample  of  the  ringing  response  would  be  perfectly  focused. 

Next  consider  the  indexing.  The  indexing  is  performed  such  that  f(i,k)  al¬ 
ways  represents  the  leading  edge  of  the  response  from  a  target  located  at 
( ifk )  regardless  of  whether  it  is  perfectly  focused  or  not.  One  ef(i,k)  is  found, 
we  now  wish  to  find  the  discrete  samples  of  the  target  ringing.  The  samples 
will  be  counted  as  m  -  0  for  the  first  sample,  m-  1  for  the  second,  and  so 
on.  We  obtain  these  samples  by  simply  incrementing  the  i  index.  The  mth 
value  is  just  /(t  +  m.k),  which  follows  from  equations  (5)  and  (7)  as 

fi,k(ma)  ~  fi-n+m,k(na)  “  f(l  +  m’k)  >  °r 

y^sj(Ti,kj  +  ma )  *  /,-„+«.*(»«)  -  /O'  +  m,k)  . 

J 

Clearly  equation  (8)  is  identical  to  equation  (1)  (i.e.s  perfect)  when  m  -  n. 
So  equation  (8)  gives  perfect  focus  at  m~n.  The  name  “short  impulse  re¬ 
sponse  approximation”  is  used  because  the  approximation  needs  to  remain 
accurate  only  over  the  duration  of  a  target’s  impulse  response. 


4.2  Examples  of  Applying  Approximation 

To  illustrate  the  use  of  equation  (8),  we  follow  these  steps: 

1.  Fix  n. 

2.  Place  a  target  (for  the  purposes  of  the  illustration)  at  say  i  =  237  in  range  on 
the  kth  bearing. 

3.  Use  equation  (7)  to  calculate  /(/, k)  for  all  (i.k). 

The  cases  of  interest  are  where  n~  0  and  where  n*  0.  Let  us  consider  each 
separately. 

4.2.1  Case  where  n  =  0 

Where  n  =  0,  use  equation  (8)  to  find  the  ringing  response  of  the  target.  The 
response  for  the  first  three  points  of  the  focused  impulse  response  is 

/237.i(0)“/23^+3,,,(0)»/z3M(0)  =  /(237^)  (case for  m-0)  , 
fi37,i(a)  -  fm-n  -  fix. TO)  -  /(238,  A)  (case  for  m  -  1)  , 

T2«)  =  /237-n  TO)  =  /237,  TO)  -  f(239,  k)  (case  for  m  -  2)  . 

Since  n  -  0,  the  amplitude  of  the  leading  edge  (m  =  0)  of  the  impulse  re¬ 
sponse  of  that  target  is  perfectly  focused.  As  m  increases,  the  approximation 
gets  worse.  So  the  approximation  is  useful  as  long  as  the  target  resonance 
dies  before  the  approximation  gets  too  bad. 

4.2.2  Case  where  n  &0 

Where  n  *  0,  perfect  focus  is  na  seconds  past  the  leading  edge  of  the  target 
impulse  response.  Suppose,  for  this  example,  n  =  2.  The  first  four  data 
points  for  the  impulse  response  are 

/237.TO)  ~/237-»+m,T«a)  -  /23,*(2g)  -  /(237,  A')  (case  for  m  -  0)  , 

/237,  T°0  “  fi3i~n  +m,  T«a)  -  /236,  T2cx)  =  Z(238-  (case  for  m  -  1)  , 

T37,  T2a)  *  fiyi-n  4 «,  T««)  -  /237.  T2a)  =  /(239,  -fe)  (case  for  m  =  2)  , 

T37.  T3a)  -  fi'si-a  4™.  T««)  -  /2JS,  T2«)  =  /(240,  *)  (case  for  m  -  3)  . 

Note  that  the  leading  edge  is  not  perfectly  focused,  as  it  was  when  n  =  0. 
The  important  point  here  is  that  one  can  choose  n  *  0  to  allow  the  leading 
edge  to  defocus  slightly  for  the  sake  of  keeping  later  points  in  better  focus. 
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5.  Error  of  Short  Impulse  Response  Approximation 


Landt,  Miller,  and  Van  Blaricum  [7]  show  the  transient  echo  response  of  a 
thin  (high-0  dipole.  Its  ringing  is  damped  after  five  cycles.  This  response 
allows  one  to  gain  an  intuitive  idea  of  the  kind  of  range  needed  in  the  ap¬ 
proximation.  A  1/2-ft  dipole  should  ring  for  five  cycles  at  1  GHz.  If  the  A/D 
sampler  collects  data  at  2  Gs/s  and  r  =  0,  then  at  m  ~  10,  data  for  all  five 
cycles  will  have  been  collected.  If  the  dipole  were  5  ft  long,  then  it  would 
ring  for  five  cycles  at  100  MHz.  So  at  m  =  100,  data  for  the  five  cycles  will 
be  collected.  Generally,  high  frequencies  damp  quickly  and  low  frequencies 
damp  slowly. 

Figure  3  illustrates  the  geometry  of  the  approximation.  The  discussion  as¬ 
sumes  this  geometry.  If  a  target  is  located  at  (i,k)  =  (237,0),  then  R  is  the 
distance  from  the  center  of  the  array  to  the  target;  a  x  is  the  distance  from  the 
end  of  the  array  to  the  target;  and  /k37  0(0)  is  the  perfectly  focused  data 
point  for  the  impulse  response  of  the  target.  The  next  point  ( m  =  1)  in  the 
impulse  response  is  approximated  by  a2  —  ax  +  ad,  or  in  general  a7  a}  + 
mad.  If  e  is  the  difference  (error)  between  the  approximation  and  the  actual, 
then 


e 


/ 


[ a ,  +  ma. 


I  T 

=  JR2  +  —  -  LR  cos  9  +  ma , 
V  2 


■  ^ j(R  +  mad)2  +  —  -  LS(R  +  mad)cos6 


(9) 


As  an  example  of  finding  the  approximation  error,  suppose  R  -  2  km,  6  = 
90°,  Ls  ~  1  km,  and  m  =  10.  How  many  degrees  off  (round  trip)  would  the 
end  points  be  at  1  GHz?  Plugging  these  numbers  into  equation  (9),  we  find 
s  =  0.0447  m.  So  the  round  trip  phase  error  at  1  GHz  is  107°,  Figure  4  is  a 
plot  of  the  error  as  a  function  of  the  position  in  the  aperture  for  0  =  90°,  76°, 
50°,  and  30°.  A  peculiar  characteristic  is  that  the  error  peaks  at  76°.  This 
peaking  is  a  result  of  working  at  a  range  of  only  twice  the  aperture  length. 
Figures  5  and  6  show  the  error  as  a  function  of  the  beam  angle  6,  for  m  =  5, 
10,  15,  and  20  at  two  ranges. 

Figure  7  is  a  plot  of  the  error  at  the  right  end  point  of  the  aperture  as  a  func¬ 
tion  of  m,  with  T—  0.  Since  the  error  at  the  right  end  is  greater  than  the  error 
at  the  left  end  for  ©angles  below  90o>  this  is  the  worst  case.  Figure  8  is  the 
same  plot  but  with  r  =  na  and  n-  4.  In  this  case,  the  leading  edge  of  the 
impulse  response  (at  m  =  0)  is  out  of  focus  by  about  45°  at  1  GHz.  Perfect 
focus  occurs  when  m  ~  n  -  4. 

Figures  4  through  8  give  an  indication  of  the  bounds  over  which  a  five- 
cycle  ring  of  a  high-0  scatterer  is  adequately  captured  by  the  approximation 
used  in  equation  (4).  These  figures  also  show  that  making  r  *  0  can  double 
the  depth  of  focus  on  resonant  targets.  The  optimum  r,  however,  depends 
on  the  frequency  and  Q  of  the  resonance. 
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1  GHz 


f>*  Efficient  Calculation 


Fast  focusing  can  be  broken  into  three  stages:  (1)  preprocessing  raw  data, 
where  interpolation  is  performed;  (2)  computing  a  set  of  polynomial  coeffi¬ 
cients  that  wall  be  used  for  fast  index  calculation;  and  (3)  performing  the  fo¬ 
cusing  summation  of  equation  (7)  using  the  coefficients  found  in  stage  2  to 
find  the  index  corresponding  to  the  proper  time  shift, 

6,1  Preprocessing 

A  method  to  efficiently  implement  time  shifting  is  needed.  The  SAR  data 
are  collected  by  an  A/D  converter  that  outputs  a  vector  of  N  numbers  (volt¬ 
ages)  at  each  aperture  position.  We  call  this  vector  i.(i').  A  time  shift  is, 
therefore,  simply  a  shift  in  the  index  i. 

Typically,  the  time  shifting  obtained  by  indexing  on  the  original  Appoint 
vector  is  not  fine  enough.  Finer  resolution  is  gained  by  a  two-stage  interpo¬ 
lator.  First,  a  high-quality  interpolator  is  used  to  produce  a  new  /fopoint 
vector  It  is  usually  implemented  by  insertion  of  Z-  1  zeros  between 
each  data  point  in  the  original  vector,  after  which  the  new  sequence  is 
passed  through  a  low-pass  FIR  (finite  impulse  response)  filter.  The  process 
results  in  a  new  vector  with  nearly  the  desired  time-shift  resolution.  The 
length,  in  this  case,  is  K  =  ZN.  We  can  gain  extra  fine  resolution  by  using  a 
floating-point  index  with  simple  linear  interpolation  to  find  a  value  between 
any  two  data  points  in  the  interpolated  data  vector.  For  example,  if  the  in¬ 
dex  i  were  6.3,  then  the  interpolated  value  would  just  be  0.7s, (6)  +  0. 3sy(7). 

The  focusing  algorithm  that  follows  uses  these  techniques  in  the  following 
sequence,  which  is  typically  referred  to  as  back  projection. 

1.  Read  in  one  A-point  vector s. (i)  of  raw'  data. 

2.  Do  an  Z-point  interpolation,  to  form  a  new  K- point  vector  Sj(i). 

3.  iterate  for  ail  pixels: 

a.  Find  the  floating  point  index  needed  for  a  pixel. 

b.  Find  the  data  value  for  that  pixel  either  by  rounding  the  index  and  grab¬ 
bing  the  value,  or  by  linearly  interpolating  between  adjacent  values. 

c.  Sum  into  that  pixel  the  data  value  obtained. 
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Fast  Index  Calculation 


6.2.1  Approach 

Typically,  the  greatest  computational  load  in  back-projection  focusing  is  in 
finding  the  index.  Let  P(i,kJ)  be  the  exact  index  needed.  Then  equation  (7), 
the  2-D  focusing  equation,  becomes 

fu‘2s^p(i’kJ))  ■  (I0) 

j 

Calculating  P(i,k,j )  exactly  involves  3-D  trigonometric  solutions  for  every 
pixel  in  the  image  at  every  position  in  the  aperture.  Such  a  calculation 
would  be  prohibitively  large.  For  practical  implementation,  a  secondary  ap¬ 
proximation  is  used  to  speed  the  computations.  Extremely  efficient  compu¬ 
tational  methods  exist  for  finding  evenly  spaced  solutions  to  polynomials. 
Therefore,  our  approach  is  to  define  a  polynomial  in  three  variables  G(i,/c,j), 
where  G(i,k,j)  ~P(i,k,f),  and  use  G  to  compute  the  index. 

The  first  issue  that  arises  is  choosing  the  order  of  the  polynomial  needed  for 
each  variable.  To  understand  the  method,  suppose  that  a  second-degree 
polynomial  is  adequate  for  each  of  the  three  variables  i, kj.  Now  the  prob¬ 
lem  can  be  restated  as  follows:  for  a  given  image  pixel  (i,k),  and  a  given 
aperture  position  (j),  find  the  coefficients  am  for  m  =  0...26  such  that 


G(i,k,j)  ~  P(i,k,j) 

-  <?(),*,/ +  w  * 

-  1C0J  +  Gik  +  c2,/2]  +  icxi  +  c4tik  +  Csjk2)i  +  [c6j  +  clyjk  +  cSJk2]i2 

“[(«o  +  aj+  a2f)  +  (a3  +  aJ  +  asfL)k+  («6  +  a~3  +aJZ)k2)  pm 
4-  [(a,,  +  a,J  +  avf)  +  (an  +  al3j  +  aiAj2)k  +  (al5  +  aj  +  a„j2)k 2 ]i 
+  f(aiH  +  al9j  +  a20j2)  +  (a2,  +  anj  4-  a2:f)k  +  (a2>  +  azJ  +  a2J2)kz]i2  . 

With  the  equation  written  in  this  format,  it  is  easy  to  see  that  one  can  code 
the  index  calculation  as  loops  nested  three  deep,  with  j  as  the  outer  loop,  k 
as  the  middle  loop,  and  i  as  the  inner  loop. 


6.2.2  Solving  for  the  Coefficients 

The  coefficients  am  can  be  found  numerically.  The  approach  is  to  calculate 
the  exact  index  required  for  M  points  in  the  image  at  each  of  H  positions  in 
the  aperture,  and  do  a  least-squares  fit  to  find  the  coefficients  am  for  m  = 
G...26.  We  can  find  a  generic  solution  to  this  problem,  for  arbitrary  image 
size,  by  using  a  pseudo- inverse  operation  to  perform  the  least-squares  fit. 
When,  matrix  A  is  not  square  but  rectangular,  then  A-1  is  known  as  the 
pseudo  inverse  and  is  defined  as 

A”3  =  (A7  *  A)"’  »  Ar ,  (12) 

An  example  explains  the  method.  To  simplify  the  example,  we  find  a  solu¬ 
tion  for  only  a  single  position  in  the  aperture  at  j  =  0.  So  we  find  the  c() 
through  c g  needed  to  calculate  the  index  needed  for  s()(G(i,kt ())).  Take  M  = 
25  points  on  a  patch  to  be  focused:  five  ranges  (close  in  to  far  out:  i  -  0,  h, 
2b,  3b,  4b)  on  each  of  five  azimuths  (left  side  to  right  side:  k  =  0,  a ,  2a,  3 a, 
4a).  The  vector  B  is  the  exact  (floating-point)  calculated  index  needed  for 
those  25  points  in  the  image.  The  matrix  A  is  set:  up  so  that  the  width  of  the 
patch  is  4a  and  the  depth  of  the  patch  is  4b.  So  we  have 


G(i,kJ)  |;._0  =  c(kq+cuqI 


+  c,  (}i  +  r, , 


+  C6.»l 


*2  r  .2  v  -y 

7,0 1  ^  *-S,0*  k 
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which  is 
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The  inverse  of  A  is  found  to  be 
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Now  C  is  calculated  as  C  =  A"'B,  and  is  used  for  finding  a  floating-point  in¬ 
dex  pointing  into  the  data  array.  This  floating-point  index  is  either  rounded, 
so  that  the  nearest  point  is  chosen,  or  is  used  for  a  linear  interpolation  be¬ 
tween  the  two  closest  data  points,  so  that  a  data  value  is  found. 

Two  options  are  available  for  finding  the  coefficients  o0  ...  a26.  The  first 
method  is  to  go  through  the  same  process  as  shown  but  for  the  full  problem. 
That  process  would  involve  enlarging  C  to  hold  the  27  coefficients  (c0  26)> 
enlarging  B  to  hold  the  ideal  values  for  more  than  one  position  in  the  aper¬ 
ture  (say  five  positions);  and  enlarging  A  to  match. 

The  second  method  is  to  solve  for  vector  C  (as  shown)  at  a  number  of  posi¬ 
tions  in  the  aperture.  Then  construct  a  polynomial  in  j  for  each  coefficient. 
This  solution  does  not  involve  inverting  the  large  A  matrix  but  will  result  in 
a  less  than  optimum  solution  in  the  least-squares  sense. 

One  advantage  of  this  second  method  is  that  it  lends  itself  to  correcting  for 
airplane  motion.  Since  A  is  independent  of  the  airplane  position,  it  is  in¬ 
verted  once.  The  matrix  multiplication  to  find  the  index  polynomial  coeffi¬ 
cients  can  be  applied  at  every  aperture  position  or  short  subaperture  if  de¬ 
sired.  The  algorithm  thus  allows  real-time  UWB  motion  compensation, 

6.2.3  Fast  Polynomial  Calculation 

Now  that  we  have  a  polynomial  to  find  the  index,  we  need  a  fast  way  to 
compute  the  polynomial.  Directly  calculating  an  Ath  degree  polynomial 
usually  requires  N  adds  and  A  multiplies.  If,  however,  a  sequence  of  solu¬ 
tions  is  desired  with  the  variable  changed  in  fixed  increments  (the  case  we 
have  here),  more  efficient  means  are  available.  An  algorithm  by  Nuttall  [8] 
describes  a  method  that  can  be  applied  here,  to  calculate  the  index  recur¬ 
sively  without  the  multiplications.  The  procedure  is  as  follows: 

1.  Let  X(f)  =  q0  +  q}i  +  q2fi  be  the  polynomial  to  be  solved,  where  i  -  0,  1,  2, 
3.... 

2.  Let  Ajfz)  =  X(i)  -  X(i  -  1)  =  -  q2+2q2i.  Observe  that  Xx{i)  ~  X 3 i  -  1)  = 

lch- 

3.  Therefore  a  recursion  can  be  set  up  where 

X, (i)  «  Xx{i  - 1)  +  2q2,  and 
X(i)-Xx(i-l)  +  Xx(i) 

The  starting  values  for  the  recursion  are 

X(0)  =  q()  and 

(0)  =  qx-q2+  2 q2i  =  qt  -  q2. 
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The  same  derivation  can  be  applied  to  an  arbitrary  degree  polynomial  to 
produce  a  recursive  formula  that  requires  N  adds  per  step.  The  technique 
can  also  be  expanded  to  multiple  dimensions  by  nesting.  Appendix  B  lists 
all  routines  used  in  the  fast  focusing  algorithm.  Subroutine  poly2  in  section 
B-5  of  appendix  B  uses  this  technique  directly. 


6.2.4  Coefficient  Generator  Program 

Figure  9  shows  how  the  image  is  broken  down,  along  with  the  names  used 
in  the  subroutines. 

The  basic  flow  chart  is  shown  in  figure  10.  In  figure  10,  the  Z  ' s  are  the  aQ 
to  a coefficients.  Each  time  the  inner  loop  goes  through  ail  the  n' s,  the 
Zp’ s  are  calculated  for  the  box  (defined  by  h  and  m)  and  the  subaperture 


Figure  9.  Image 
partitioning  and 
notation. 


k  boxes  =  3 


Figure  10.  Coefficient 
generator  flow  chart. 
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(defined  by  I).  The  flow  chart  shows  that  the  polynomial  degree  for  the  Cn 
is  not  always  the  same.  The  degree  was  optimized  based  on  the  errors 
caused  when  a  lower  degree  was  tried.  Appendix  A  provides  a  complete 
listing  of  the  program  that  generates  the  focusing  coefficients  from  an  input 
file  with  the  geometry. 

Fast  Focusing  Algorithm 

The  fast  focusing  algorithm  simply  applies  the  fast  polynomial  solution 
technique  to  find  the  indexing  and  does  the  signal  summation  into  all  the 

is  done  in 
jf subrou- 
coded  ex¬ 
technique 
3r  illustra- 
ken  into  a 
fficients. 
derivation 
us  a  sinele 


A8.1  jEIISij  ssjfp1 


liliTj  JgiTil'l 


Inner_Loop  (r,  s,  s_max,  q 


float  *q; 


int  s_max; 


int  i_pixes; 


t__stop_i  =  i  +  j_pixes-l; 

f-q[2]+q[2]; 

Rl=q[l]-q[2j; 

RO  =  q[0]; 


if  RU  <  -0.5  then  repeat 
{f++; 

R1  =  Rl+T; 

RO  =  R0+R1 }; 
until  RO  >-0.5}; 

If  (f  >  fjstogj)  then  return; 


index  ~  Round(RO); 
*f  =  *f  +  *(s+index); 


repeat 


mdex  =  Roundfl 
if  index>s_max 


until  (f  =  f_siop); 
return; 


Pomter  to  first  pixel  in  a  line;  f(i  start, k) 


vector 


Pointer  to  coefhc  tents  vector  for  index  calculation 


s  range  is 


;  along  i  axis  is  i  nixes 


up  initial  co 


Perform  clipping  when  index  <  0  (before 
data). 


4~4“*  j 

■J 

ncrement  pointe 

r  to  next 

pixe 

i. 

UR 

1+T; 

m 

ncrement  index  RO  = 

R0+R1. 

when  index  is  beyond  actual  data  record  length. 


Sum  into  pixel  the  new  data 


Nesting  of  the  recursive  approach  means  that  the  coefficients  q0,  q{ ,  and  q2 
are  each  formed  by  polynomials  in  k.  Each  polynomial  is  incremented  re¬ 
cursively  in  the  subroutine  MiddleJLoop  (display  2).  Thus,  MiddleJLoop 
sums  into  a  single  box  the  data  from  a  single  aperture  position. 


Display  2.  Subroutine  MiddleJLoop 
"Code" 


MiddIe_Loop{f,  s,  e,  sjnax, 

ijjixesTkJnc,  box_offset); 

___  _______ 


float  *s; 
float  *c; 


Comments 


Middle  loop  subroutine 


Pointer  to  first  pixel  of  box 


int  s_max; 
int  kjnc; 
int  ijpixes; 


int  box_pffset; 
f  last  -  f  +  box  offset; 


qO  T  =  c[2]+ cf21; 
qO_Rl  *  c[l  J — c[2j; 
q[0]  =  c[G]; 


ql  T  =  cT5]+c{5]; 
qlJRl  =  c(4]-c[5]; 


q[il 


q2  T  =  c[S]+c[8j; 
q2_Rl  =  c[7]-cf8]; 

q[2]  ^  c[6j; _ _ 

Call  Inner_Loop(f,s,q..); 
repeat 


f  =  f  +  k  inc; 


Pointer  to  jih  data  vector 


Pointer  to  coefficients  vector  for  index  calculation 
s  range  Is  s[0]...s[s_max] 


F[i,k]  =  *(f  +  i  +  k  *  kjnc)  so  k_inc  is  total  range  line  length 


Size  of  the  box  (pixels)  along  i  axis. 


Number  to  add  to  /  to  move  pointer  to  last  line  in  box. 


fjast  =  address  of  the  first  pixel  in  last  line  of  box. 
Set  up  initial  conditions  for  q<> 


Set  up  initial  conditions  for  qx. 


Set  up  initial  conditions  for  q2. 


Perform  summation  on  first  line. 
Start  Iood  for  rest  of  lines. 


Increment  pointer  to  next  line  in  box. 


qO_Rl  «qO_Rl+qO_T; 
q[0]  =  q[0]+qOJRl; 


Iterate  q{> 


qlJRl  =ql_Rl+ql_T; 


Iterate  qh 


q2_RI  =  q2_Rl+q2Jf; 
q[2]  =  q[2]+q2JRl; 


Iterate  q2. 


Cal)  InnerJLoop(..); 


Perform  summation  on  current  line  in  box. 


Display  3.  Subroutine  Outer_Loop  (cont’d) 

Code 

Comments 

END  (for  ijbox); 

!  Loop  through  all  boxes. 

END  (for  kjx>x); 

repeat 

Loop  through  all  aperture  positions. 

Increment  aperture. 

read  s[0..sjmax]; 

Position  read  data  for  that  aperture  position 

f  ss  fo  —  (box_inc_k+box_inc_i); 

Initialize  /. 

box  =  -2; 

for  k_box  =  0  to  k_boxes-l; 
box  =  box  +  1; 
f  =  f  +  box_inc_k; 
for  i_box  =  0  to  i_boxes-l; 
box  =  box  +  1; 
f  =  f  +  box_JncJ; 

Set  up  to  loop  through  all  boxes. 

cO_Rl[box]  =  cO_Rl[box]+cO_T[box]; 
c[0,box  j  =  c[0,box]+c0_R1  [box]: 

Iterate  c0  for  current  box. 

cl_Rl[box]  =  cl_Rl[box]+cj_T[boxj; 
cjl.box]  =  cfl,bojt]+cl_Rl[box]; 

Iterate  C\  for  current  box. 

Iterate  for  current  box. 

c8_Rl[boxj  =  c8_Rl  [box]+c#_T[box]; 
c[8,box]  =  c[8,box]+c8JRl[box]; 

Iterate  eg  for  current  box. 

Cali  M iddle_Loop{f,s,c(box],. .}; 

Sum  into  current  box. 

END  (for  i_box); 

END  (for  k_box); 

Loop  until  all  boxes  done. 

Until  j  =  j_stop; 

Return  (END  Outer_Loop) 

Loop  until  all  aperture  positions  are  done. 

6.4  Considerations  for  Fast  Focusing  Polynomial  Order  Selection 

A  study  was  conducted  to  determine  the  effect  of  polynomial  order  on  the 
image  size  and  aperture  length.  The  worst-case  position  in  the  aperture  (j  = 
3  in  fig.  3,  p  17)  and  the  worst-case  area  in  the  image  ((i.k)  =  (237,3)  in  fig. 
3)  were  chosen  for  the  analysis.  The  sample  rate  was  4  Gs/s  with  a  record 
length  (s_max)  of  4096  points.  The  elevation  was  60  ft;  otherwise  the  ge¬ 
ometry  was  as  in  figure  3,  with  an  aperture  2,  ?  =  384  ft,  squint  angle  0*  85°, 
and  slant  range  R  -  550  ft.  If  the  error  in  the  floating-point  index  is  re¬ 
stricted  to  ±0.5,  then  this  geometry  produces  the  following  results: 


Polynomial 

Maximum 

Maximum 

degree 

i_j>ixes 

kjjtxes 

1 

228 

27 

2 

964 

105 

% 

2191 

315 

4 

— 

657 

By  using  the  nested  recursive  technique  as  shown  in  the  example  subrou¬ 
tines,  we  compute  the  computational  load  as  shown  in  table  2. 


to  t> 


urouping  terms,  we  get 

L  =  (inner  loop  adds)  +  (middle  loop  adds)  -f  (outer  loop  adds) 

=  IKJB(i'  +  1)  +  KJBii'  +  1  )k'  +  JB(i'  +  1  )(k'  +  1);' . 

From  this  equation,  it  is  clear  that  a  low-order  polynomial  is  crucial  in  th 
inner  loop.  Conversely,  the  order  of  the  polynomial  in  the  outer  loop  can  b 
high  with  little  impact  on  the  overall  speed.  Table  3  shows  the  eomputa 
tional  load  for  various  configurations  with  j'  =  3  and  J  =  2304. 


Table  2.  Computational 
load. 


Symbol  Definition  _ _ 

i  Polynomial  order  for  i  index 
k'  Polynomial  order  for  k  index 
j'  Polynomial  order  for  j  index 
J  Number  of  points  in  the  aperture 

1  Number  of  range  bins  (pixels)  in  box  (i_pixes) 

K  Number  of  azimuth  bins  (pixels)  in  box  (k_pixes) 

5  Number  of  boxes 

L0  l(i'  +  1);  adds  per  Inner_Loop  call 

Lj  K[Lq  +  (C  +  1)A(J;  adds  per  Middle_Loop  call 

L  J[Lj  +  (/'  +  !)(&'  +  I);"];  adds  per  Outer_Loop  (adds  per  image) 


ill 


Table  3.  Matrix  showing  computational  load  versus  partitioning  for  j'  -  3  and  J  =  2304. 


k' 

1 

2 

3 

4 

K 

27 

103 

256 

512 

\ 

Boxes  across 

19 

5 

2 

1 

r 

/ 

Boxes 

down 

\jbtai  pixels 
N\  across 
ToiaF’v 
pixels  \ 
down  \ 

513 

515 

512 

512 

i 

227 

18 

4086 

1  =  9.711 

S  =  21,888 
B  -  342 

L  =  9.786 
S  =  8640 

5  =  90 

£  =  9.769 
5  =  4608 

5  =  36 

£  =  9.810 

5  =  2880 
5=  18 

2 

817 

5 

4085 

L  =  14.51 

5*  =9120 

5  =  95 

L  -  14.58 
5  =  3600 

5  =  25 

L  =  14.51 
5=  1920 
5=  10 

L  =  14.53 

5  =  1200 
5=5 

3  2043 

2 

4086 

L  -  19.33 

5  =  4864 

5  =  38 

£.=  19.41 
5=  1920 

5  =  10 

£=19.31 

5  =  1024 

5  =  4 

£  =  19.32 

5  =  640 

5  =  2 

=  total  number  of  adds  in  giga-adds. 

=  total  number  of  boxes  in  the  image. 

S  =  bytes  of  storage  needed  for  the  coefficient  tables  ~  (i '  +  l)(k'  +  l)(j'  +  1}B(4  bytes  Sword). 
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Implementation  Notes 

The  algorithm  developed  here  requires  no  multiplication  except  in  the  pre¬ 
processing  (interpolation)  stage.  There  already  exist  DSP  (digital  signal 
processor)  chips  whose  architecture  is  ideal  for  the  FIR  interpolation  task. 
The  bulk  of  the  computations,  however,  occur  in  the  summing  and  index 
calculation  stage.  Two  facts  make  the  summing  algorithm  attractive.  First, 
adders  take  considerably  less  area  to  implement  in  a  VLSI  (very-large-scale 
integration)  chip  than  multipliers.  Second,  parallel  operation  is  extremely 
simple;  the  image  can  simply  be  broken  into  boxes  with  a  separate  proces¬ 
sor  working  independently  on  each  box.  In  this  case,  only  the  coefficient 
table  would  be  different  between  processors,  and  a  broadcast  mode  would 
be  needed  to  pass  the  signal  data  to  each  processor.  Therefore,  design  of  a 
custom  LSI  chip  appears  practical  and  could  result  in  real-time  speeds. 

If  an  off-the-shelf  DSP  chip  with  a  parallel  adder  and  multiplier  is  used  to 
perform  this  algorithm,  then  other  options  become  available.  For  example, 
the  multiplier  can  be  used  in  the  summing  algorithm  to  do  interpolation  on 
the  indexing  rather  than  rounding.  The  multiplication  required  can  be  done 
during  other  add  cycles  so  that  it  takes  zero  time.  Another  possibility  is  to 
calculate  the  index  polynomial  directly  instead  of  recursively.  Finally,  since 
all  the  operations  can  be  fixed  point,  the  address-generator  ALU  (arithmetic 
logic  unit)  can  sometimes  be  used  in  parallel  with  the  floating-point  units. 


7.  Beam  Patterns  and  Sidelobe  Structure 


To  develop  an  intuitive  understanding  for  the  beam  shape,  we  assembled  16 
figures  to  display,  in  the  time  domain,  the  main-lobe  and  sidelobe  beam  pat¬ 
terns.  We  simulated  a  resonant  target  by  generating  a  data  record  for  every 
position  in  the  aperture.  The  simulated  scene  was  a  single-point  target.  The 
target  bearing  (see  fig.  3)  was  squinted  15°  off  broadside  (6  =  105°).  Range 
to  the  target  Rwas  750  ft,  the  aperture  length  Ls  was  385  ft,  and  the  aperture 
height  was  60  ft.  The  point  target  had  an  impulse  response  of 


s(t)> 


sin(2;i//)[0.5  +  0.5cos(4jrf/)]  for0<f< 


4  f 


sin(2jr/i)exp 


4/ 


0.23/ 


for  t  & 


j_ 

4/ 


(15) 


The  data  simulated  a  2-GHz  sample  rate.  A  record  length  of  2048  samples 
was  made  for  each  aperture  position.  Each  record  was  preprocessed  with  a 
times-8  interpolator,  producing  16K  records.  Interpolation  was  done  by  the 
standard  FIR  filter  method.  A  255-tap  Parks-McClellan  low-pass  filter  with 
a  950-MHz  cutoff  was  used  to  do  the  interpolation.  Equation  (7)  (with  n  = 
0)  was  used  to  produce  the  focused  beams.  A  Hilbert  transform  was  used  to 
obtain  the  magnitude  that  is  plotted  in  the  figures.  Plots  were  made  on  a 
decibel  scale.  Four  3-D  plots  were  made,  showing  two  viewing  angles  (on 
axis  and  above  axis)  and  two  aperture  weighting  functions  (Hamming  and 
rectangular)  for  each  of  four  targets  with  different  resonant  frequencies  (50, 
200,  400,  and  900  MHz).  These  are  shown  in  figures  11  to  26. 

The  on-axis  plots  were  made  so  that  amplitude  values  could  be  easily  read 
off  the  plots.  The  above-axis  plots  were  made  to  reveal  the  sidelobe  struc¬ 
ture  and  the  ring-down  time  of  the  target.  The  axis  labels  were  shifted  so 
that  0  range  was  at  the  point  target  and  0°  was  the  bearing  centered  on  the 
target.  Notice  that  all  plots  are  3-D,  even  the  on-axis  plots.  The  various 
horizontal  lines  in  the  on-axis  plots  are  the  beam  pattern  on  the  Ith  range 
bin.  These  horizontal  lines  are  identical  to  those  shown  in  the  above-axis 
plots;  they  are  just  being  viewed  “end  on.” 

These  figures  are  unique  in  that  they  show  how  the  sidelobes  spread  in  time 
as  one  moves  off  the  main  beam.  Because  of  the  sharp  rise  time  of  the  tar¬ 
get  echo,  the  contribution  from  the  near  and  far  ends  of  the  aperture  can  be 
seen  as  the  early  and  late  peaks  in  the  sidelobe  structure.  Hamming  weight¬ 
ing  was  applied  to  the  array  so  that  one  can  see  the  effect  of  weighting  on 
the  sidelobe  structure.  It  is  interesting  to  note  that  although  the  peak 
sidelobe  levels  drop  only  marginally,  the  average  or  integrated  sidelobe  lev¬ 
els  are  significantly  lower  when  tapered  aperture  weighting  is  used.  The 
weighting  also  reduces  the  near  and  far  peaks  in  the  sidelobes  since  the 
weighting  reduces  the  contribution  of  the  array  ends. 
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Since  the  antenna  is  a  linear  system,  it  is  no  surprise  that  the  beam  pattern 
behaves  basically  as  classic  antenna  theory  predicts.  The  aperture  beam  pat¬ 
tern  is  a  function  of  (1)  how  long  the  aperture  is  relative  to  wavelength,  and 
(2)  what  kind  of  weighting  is  used  to  sum  the  points  in  the  aperture.  If  uni¬ 
form  aperture  weighting  is  used  (a  straight  summation),  then  the  typical 
sin (ftip)l(fi\p)  pattern  occurs,  where  ip  is  the  angle  off  the  main  beam,  and  fi 
cc  L/k.  f3  establishes  the  width  of  the  main  beam.  As  the  wavelength  gets 
small,  or  as  the  aperture  gets  long,  the  beam  gets  narrow. 


Figure  11.  50-MHz 
point  reflector,  equal 
weighting,  end  view. 
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Figure  15.  200-MHz 
point  reflector,  equal 
weighting,  end  view. 
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Figure  17. 200-MHz 
point  reflector,  equal 
weighting,  above  axis, 


Figure  18. 200-MHz 
point  reflector, 
Hamming  weighting, 
above  axis. 


Figure  19.  400-MHz 
point  reflector, 
Hamming  weighting, 
end  view. 
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Figure  20.  400-MHz 
point  reflector, 
Hamming  weighting, 
end  view. 
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Figure  23.  900-MHz 
point  reflector,  equal 
weighting,  end  view. 
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Figure  24.  900-MHz 
point  reflector, 
Hamming  weighting, 
end  view. 
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Figure  25.  900-MHz 
point  reflector,  equal 
weighting,  above  axis. 


Figure  26.  900-MHz 
point  reflector, 
Hamming  weighting, 
above  axis. 


Image  Amplitude  fdB)  Image  Amplitude  (dB) 


8.  Summary 


This  report  identifies  a  frequency-independent  processing  algorithm  to  get  a 
perfectly  focused  impulse  response  from  any  object  at  any  position  from  an 
ultra-wide-bandwidth  synthetic  aperture  radar.  The  report  also  presents  an 
approximation  that  reduces  the  computational  complexity  of  the  algorithm. 
An  error  analysis  of  this  approximation  demonstrates  its  applicability  to  fo¬ 
cus  even  high-g  objects  in  the  near  field  of  an  aperture.  An  implementation 
that  is  both  computationally  efficient  and  applicable  to  real-time  motion 
compensation  is  also  presented.  Finally,  plots  are  presented  demonstrating 
the  capability  of  the  algorithm  to  focus  ringing  targets  over  an  18-to-l 
bandwidth. 
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Appendix  A 


A-l.  Main  Program  to  Calculate  Coefficients  (main.c) 

This  is  a  listing  of  the  code  that  calculates  the  coefficients  needed  by  the 
fast  focusing  algorithm.  It  takes  as  inputs  the  geometry  of  the  synthetic  ap¬ 
erture  radar  (SAR)  and  image  area,  and  it  outputs  the  coefficient  vectors 
required. 


Program  :  coefgen 

Description  :  This  is  the  main  program  to  generate  coefficients 
for  fast  focusing  algorithm  . 


a**************************************************************************** 


#include  <malloc.h> 

#include  <stdio.h> 

#include  <math.h> 

#include  "coefgen. const" 

#include  "coefgen.var" 

main() 

{ 

long  cacheset; 

FILE  *fp,*tstart_fp; 
char  fn[80]; 
long  i; 

ia=(BOX__SIZE_AZIMUTH-4)/4; 

ib=(BOX_SIZE_RADIAL-3)/4; 

/*  initialize  variables  */ 
coefgen_varinit(); 

/*  initialize  pseudo  inver  matrix  */ 
initmatO; 

printf("Enter  coefficient  a  output  file  name  ===>  "); 

scanf("%s",fn); 

fp=fopen(fn,  "w"); 

if(fp==NULL) 

{ 

printf("Error  open  output  file  \n"); 
exit(l); 

} 

/* 
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printf("Enter  data  acquisition  timing  file  name  ===>  "); 
scanf("%s",fn); 

7 

/* 

tstart_fp=fopen("deIaytime.dat'',Y'); 

if(tstart_fp==NULL) 

{ 

printf("Error  open  input  file  \n"); 
exit(O); 

} 

read_delaytime(tstart_fp); 

7 

coefgen(fp); 

} 


A-2.  Main  Coefficient  Generator  Routine  (coefgen.c) 

#include  "coefgen.h"  /^geometry  file*/ 

#include  "coefgen.var"  /*list  of  variables7 

#include  <stdio.h> 

#include  <math.h> 
extern  index(); 
extern  mxmul_Q; 
extern  poly_fit(); 

/*  Note  that  Reference  point  is  now  taken  at  coordinate  kpixel=  n_azimuth/2-l  7 

/*  ipixel=  n_radia!/2  -1  7 

/*  At  every  position  data  were  taken  so  that  the  reference  point  is  at  center  */ 

/*  data  buffer  7 

void  coefgen(fp) 

FILE  *fp;  / 

{ 

float  input[NPOINT_APER/4]; 
long  i,j; 

long  section, kbox,ibox,position,position_start,posi  tion_stop; 
long  ipoint,kpoint; 

for(i=0;i<NPOINT_APER/n_section;i++) 

{  . " 

input[i]=(float)i; 

} 

for(section=0;section<n_section;section++) 

{ 

position_start=section*n_aper/n_section; 

position_stop=position_start+n_aper/n_section; 
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for(kbox=0;kbox<NBOX_AZIMUTH;kbox++) 

{ 

for(ibox=0;ibox<NBOX_RADIAL;ibox++) 

{ 

printf("processing  section  %d  ibox  %d  kbox  %d\n", section, ibox,kbox); 
for(position=position_start;position<position_stop;position++) 

{ 

/***** *Calculating  Actual  index  for  all  samples  points  in  a  box**********/ 
i=0; 

for(ipoint=0;ipoint<NPOINT_BOXSAMPLE_RADlAL;ipoint++) 

{ 

for(kpoint=0;kpoint<NPOINT_BOXSAMPLE_AZIMUTH;kpoint++) 

{ 

ipixel=ibox*BOX_SIZE_RADIAL+ipoint*ib; 

jpixel=kbox*BOX_SIZE_AZIMUTH+kpoint*ia; 

d=((float)position-((float)n_aper/2.-l.))*(aper_length/((f!oat)n_aper-l.));  /*  distance  from  radar 
to  center  of  aperture  */ 

x=sqrt(rcenter_ref*rcenter_ref-radar_height2); 

r_ref=sqrt(radar_height2+x*x*c_theta_center_ref*c_theta_center_ref+(d- 
x  *  s_theta_ce  n  ter_re  f)  *  (d-x  *  s_t  heta_ce  n  ter_ref)); 
rrnin=r_ref-d_rcenter*((float)n_radial/2.-l.); 

index(&ipixel,&jpixel,&findex[i]); 

i++; 

} 

} 

/*  generate  coef  for  this  box  at  this  position  */ 

mxmuls(mat,&mat_c_stride,&stride,findex,&findex_c_stride,  Astride,  coef,&coef_c_stride, 
&stride,&n_coef,&i_one,&n_boxsample); 

/*  save  in  coefc[i][position]  */ 

for(i=0;i<COEF_SIZE;i++) 

{ 

coefc[i][position-position_start]=coef[ij; 

} 

/*  for  debug  only  */ 

/*  printf("pos=  %d  tstart=  %,3e\n",position,2.0*rmin/c);  */ 

/*  end  debug  */ 

}  /*  position  */ 
for(i=0;i<COEF_SIZE;i++) 

{ 

printf("Doing  polyfit  for  coef  %d\n",i); 
poly__fit(input,&coefc[i][0],n_aper/n_section,deg[i],&coefa[i][0]); 
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for(j=0;j<(deg[i]+ 1  );j++) 

{ 

fprintf(fp,"%e\n",coefa[i][ij);  /*  save  a[ij  to  file  */ 

} 

} 

}  /*  ibox  */ 

}  /*  kbox  */ 

}  /*  section  */ 

}  /*  subroutine  */ 


A-3.  Subroutine  Find  Ideal  Index  Vector  (index.c) 


Subroutine: 

index.c 

Input: 

ipixel 

pixel  radial  coordinate 

jpixel 

pixel  angle  coordinate 

a__span 

angle  spanned  by  the  patch 

a__ofset 

offset  angle  of  the  center  line 

d 

distance  from  middle  position 

d__rcenter 

sampling  range 

dr8 

(sampling  range )/8 

rmincenter 

min  range  rom  center  position 

hgt 

height  of  building 

hgt2 

square  the  height 

length 

length  of  the  building 

Output: 

findex 

floating  point  index  to  ivalue  of  that  pixel 

************************************************************************ 
#include  "coefgen.const" 

#include  "coefgen.var" 

#include  <math.h> 

index(ipixel,kpixel,findex) 
int  *ipixel,*kpixel; 
float  *findex; 

{ 


theta=a_ofset-a_span/2.+(f!oat)*kpixel*d_theta; 

c_theta=cos(theta); 

s_theta=sin(theta); 

rcenter=rmincenter+(float)*ipixel*d_rcenter; 
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x=sqrt(rcenter*rcenter-radar_height2); 

r=sqrt(radar_height2+x*x*cjheta*c_theta+(d-x*s_theta)*(d-x*s_theta)); 

*findex=(r-rmin)/drS; 

} 


A-4.  Geometry  File  —  Declare  All  Global  Constants 
(coefgen.h) 

#define  NPOINT_AZIMUTH  600  /*  number  of  bearing  lines  7 

//define  NPOINT_RADIAL  4095  /*  number  of  radial  lines  7 

#define  NPOINT_APER  2304  /*  #  of  positions  in  aperture  7 

#define  NPOINT_DATA  2048  /*  number  of  original  data  points  7 

//define  NPOINT_DATA_INTER  NPOINT_DATA*8 /*  number  of  data  points  after 

interpolated  7 

//define  PI  3.141592654 

#define  RADARJHEIGHT  60.  /*  radar  height  in  feet  7 

#define  A_OFSET  -15.0  /*  offset  angle  from  center  line  7 

#define  A_SPAN  54.0  /*  spanning  angle  of  the  patch  7 

#define  RMINCENTER  550.  /*  range  from  center  position  to 

nearest  point  on  the  patch  7 

//define  RCENTERREF  802.0  /*  range  (ft)  of  ref.  point  from  center  pos.  7 

#define  THETA_CENTER_REF  -15.0  /*  angle  (deg)of  ref.  point  from  center  pos.  7 
#define  SAMPLING_RATE  2.0e09  /*  sampling  frequency  of  signal  7 

#define  SAMPLING  PERIOD  5.0e-10  /*  ts=l/fs  7 


#define  BOX_SIZE__RADIAL  195 
#define  BOX_SIZE_AZIMUTH  100 

#define  NB O X_R ADI AL  NPOINT_RADIAL/BOX_SIZE_RADIAL 
#define  NBOX  AZIMUTH  NPOINT_AZIMUTH/BOX_SIZE_AZIMUTH 

#define  SECTION_SIZE  4  /*  #  of  sections  for  one  aperture  */ 

#define  COEF_SIZE  6  /*  #  of  coeffs  for  curve  fit  */ 

#define  MAXDEG  3  /*  maximum  degree  for  poly,  fit  */ 

#define  NPOINT_BOXSAMPLE_RADIAL  5  /*  #  of  samples  in  a  box  in  radial 

direction  7 

//define  NPOINT_BOXSAMPLE_AZIMUTH  5  /*  #  of  samples  in  azimuth  direction 

for  every  box  7 

//define  NPOINT_BOXSAMPLE 

NPOINT_BOXSAMPLE_RADIAL*NPOINT_BOXSAMPLE_AZIMUTH 
/*  #  of  samples  in  a  box  for  curve  fit  7 
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A-5.  Declare  All  Global  Variables  (coefgen.var) 

struct  complex  {  float  r,i;};  /*  define  a  complex  type  */ 

/*  Data  Variables  */ 

float  pixel [NPOINT_AZIMUTH][NPOINT_RADIAL];  /*  array  of  image  */ 

float  data[NPOINT_DATA];  /*  original  data  bufer  */ 

float  data_inter[NPOINT_DATA_INTER];  /*  interpolated  data  buffer  */ 

struct  complex  data_inter_fft[NPOINT_DATA_INTER/2];  /*  Real->Complex  Forward  FFT 

of  data_inter  */ 

float  filter_coef[NP01NT_DATAJNTER];  /*  filter  coefficients  */ 

struct  complex  filter Jft[NPOINTJ3ATA_INTER/2];  /*  Real->CompIex  Forward  FFT 

of  filter_coef  */ 

float  nyquist_filter;  /*  value  of  FFT  of  Filter  at  Nyquist  point  */ 

float  nyquist__data;  /*  value  of  FFT  of  data  at  Nyquist  Point  */ 


/*  Geometry  Variables  */ 

float  a_ofset;  /*  offset  angle  from  the  center  line  */ 

float  a_span;  /*  spanning  angle  of  the  patch  */ 

float  d_theta;  /*  delta  angle  */ 

float  x, theta;  /*  coordinate  of  a  pixel  */ 

float  c_theta,s_theta;  /*  cosine  and  sine  of  theha  */ 


float  d;  /*  distance  from  radar  to  center  position 

positive  to  the  right,  neg.  to  the  left  */ 
float  reenter, rmincenter,rmaxcenter;  /*  range  from  center  position  */ 
float  rcenter_ref;  f*  range  from  center  pos.  to  ref  point  */ 

float  theta_center_ref;  /*  angle  of  ref.  point  from  center  position  */ 


float  rjref;  /*  range  from  any  pos.  to  ref.  point  */ 

float  d_rcenter;  /*  sampling  distance  from  center  position  */ 

float  dr8;  /*  (sampling  distance)/8  */ 

float  r,rmin,rmax;  /*  range  from  any  position  */ 

float  aper_length;  /*  length  of  aperture  */ 

float  radar_height;  /*  height  of  radar  */ 

float  radar_height2;  /*  square  the  height  of  radar  */ 

/*  Radar  Variables  */ 

float  c;  /*  speed  of  wave  */ 

float  ts;  /*  sampling  period  of  signal  */ 

float  fs;  /*  sampling  frequency  of  signal  */ 

long  n_azimuth;  /*  #  of  points  in  azimuth  direction  */ 

long  n_radial;  /*  #  of  points  in  radial  direction  */ 

long  n_aper;  /*  #  of  positions  in  aperture  */ 

long  pix_size;  /*  #  of  points  in  pixel  array  */ 
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long  n_data;  /*  number  of  original  data  points  */ 

long  n_data_inter;  /*  number  of  interpolated  data  points  */ 

long  n_data_inter_half;  /*  1/2  #  of  interpolated  data  points  */ 

long  n_boxsamp!e;  /*  #  of  samples  in  a  box  for  curve  fit  */ 

long  n_coef;  /*  #  of  coefficients  curve  fit  */ 

long  n_section;  /*  #  of  sections  for  one  aperture  V 


/*  SSL  VAriables  */ 

float  f_zero;  /*  floating  point  zero  */ 

float  f_one;  /*  floating  point  1  */ 

long  i_one;  /*  integer  1  */ 

long  stride;  /*  stride  for  supercard  ssl  */ 


/*  working  variables  */ 
long  n_position; 

float  e_index,a_index, error, m axe rror,m inerror; 

1  ong  i e rr_m ax ,j err_max,ierr_m i n ,j er r__m i n ; 
long  column_max,row_max,column_min,row_min; 
long  ia,ib; 
long  ipixel,jpixel; 

float  c_theta_center_ref,s_theta_center_ref; 
float 

mat[COEF_SIZE][NPOINT_BOXSAMPLE],findex[NPO!NT_BOXS  AMPLE], coef[COEF_SIZE]; 
float  coefc[COEF_SIZE]  [NP01NT_APER/4];  /*  c  coef.  for  each  section  */ 
float  coefa[COEF_SIZE][MAXDEG+l];  /*  a  coef.  for  each  c  */ 

long  deg[COEF_SIZE];  /*  degree  for  each  coefc  */ 

long  mat_c_stride,findex_c_stride,coef_c_stride; 

float  tstart[NPOINT_APER];  /*  time  from  trigger  to  start  data  acq.  at  */ 

/*  each  position  */ 

A-6.  Matrix  Pseudo  Inverse  Coefficients  (initmat.c) 

#include  <stdio.h> 

#include  <math.h> 

#include  "coefgen. const" 

#include  "coefgen. var" 

initmatO 

{ 

int  i,j; 
float  a,b; 

a=((float)BOX_SIZE_AZIMUTH-4.)/4.; 

b=((float)BOX_SIZE_RADIAL-3.)/4.; 
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mat[0][0]=93./175.; 

mat[0][l]=27./175.; 

mat[0][2]=(-9.)/l  75.; 

mat[0][3]=(-3.)/35.; 

mat[0][4]=9./175.; 

mat[0][5]=62./175.; 

mat[0][6]=l8./175.; 

mat[0][7]=(-6.)/175.; 

mat[0][8]=(-2.)/35.; 

mat[0][9]=6./175.; 

mat[0][10]=31./175.; 

mat[0][ll]=9./175.; 

mat[0][12]=(-3.)/175.; 

mat[0][13]=(-l.)/35.; 

mat[0][14]=3./175.; 

mat[0][15]=0.; 

mat[0][16]=0.; 

mat[0][17]=0.; 

mat[0][18]=0.; 

mat[0][19]=0.; 

raat[0][20]=(-31  .)/l75.; 

mat[0][21]=(-9.)/175.; 

mat[0][22]=3./175.; 

mat[0][23]=l./35.; 

mat[0][24]=(-3.)/175.; 

mat[l][0]=(-81.)/(175.*a); 

mat[l][l]=39./(350.*a); 

mat[l][2]=l2./(35.*a); 

mat[l][3]=81./(350.*a); 

mat[l][4]=(-39.)/(175.*a); 

mat[l][5]=(-54.)/(175.*a); 

mat[l][6]=l3./(175.*a); 

mat[l][7]=8./(35.*a); 

mat[l][8]=27./(175.*a); 

mat[l][9]==(-26.)/(175.*a); 

mat[l][10]=(-27.)/(l75.*a); 

mat[l][ll]=13./(350.*a); 

mat[l][12]=4./(35.*a); 

mat[l][13]=27./(350.*a); 

mat[l][14]=(-13.)/(175.*a); 

mat[l][15]=0.; 

mat[l][16]=0.; 

mat[l][17]=0.; 

mat[l][18]=0.; 
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mat[l][19]=0.; 

mat[l][20]=27./(l75.*a); 

mat[l][21]=(-13.)/(350.*a); 

mat[l][22]=(-4.)/(35.*a); 

mat[l][23]=(-27.)/(350,*a); 

mat[l][24]=13./(175.*a); 

mat[2][0]=3./(35.*pow(a,2.)); 

mat[2][l]=(-3.)/(70.*pow(a,2.)); 

mat[2]  [  2]  =(-3 .)/ (35 .  *  pow(a,2 .)); 

mat[2][3]=(-3.)/(70.*pow(a,2.)); 

mat[2][4]=3./(35.*pow(a,2.)); 

mat[2][5]=2./(35.*pow(a,2.)); 

mat[2][6]=(-l.)/(35.*pow(a,2.)); 

mat[2][7]=(-2.)/(35.*pow(a,2.)); 

mat[2][8]=(-l.)/(35.*pow(a,2.)); 

mat[2][9]=2./(35.*pow(a,2.)); 

mat[2][10]=l  ./(35.*pow(a,2.)); 

mat[2][l  l]=(-l.)/(70.*pow(a,2.)); 

mat[2][12]=(-L)/(35.*pow(a,2.)); 

mat[2][13]=(-l.)/(70.*pow(a,2.)); 

mat[2][14]=l./(35.*pow(a,2.»; 

mat[2][15]=0.; 

mat[2][16]=0.; 

mat[2][17]=0.; 

mat[2][18]=0.; 

mat[2][19]=0.; 

mat[2][20]=(- 1  .)/(35.*pow(a,2.)); 
ma  t[2][2 1  ]= 1  ./(70.  *pow(a,2.)); 
mat[2][22]=l./(35.*pow(a,2.)); 
mat[2][23]=l./(70.*pow(a,2.)); 
mat[2][24]-(- 1  ,)/(35.  *pow(a,2.)); 

mat[3][0]=(-31.)/(l'75.*b); 

mat[3][l]=(-9.)/(l75.*b); 

mat[3][2]=3./(175.*b); 

mat[3][3]=l./(35.*b); 

mat[3}[4]=(-3.)/(l75.*b); 

mat[3][5]=(-31.)/(350.*b); 

mat[3][6]=(-9.)/(350.*b); 

mat[3][7]=3./(350.*b); 

mat[3][8]=l./(70.*b); 

mat[3][9]=(-3.)/(350.*b); 

mat[3][10]=0.; 

mat[3][ll]=0.; 
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mat[3][l2]=0.; 

mat[3][13]=0.; 

mat[3][14]=0.; 

raat[3][l5]=31./(350.*b); 

mat[3][16]=9./(350.*b); 

mat[3][17]=(-3.)/(350.*b); 

mat[3][l8]=(-l.)/(70.*b); 

mat[3][19]=3./(350.*b); 

mat[3][20]=31./(175.*b); 

mat[3][21]=9./(175.*b); 

mat[3][22]=(-3.)/(175.*b); 

mat[3][23]=(-l.)/(35.*b); 

mat[3][24]=3./(175.*b); 

mat[4][0]=27./(175.*a*b); 

mat[4][l]=(-13.)/(350.*a*b); 

mat[4][2]=(-4.)/(35.*a*b); 

mat[4][3]=(-27.)/(350.*a*b); 

mat[4][4]=13./(175.*a*b); 

mat[4][5]=27./(350.*a*b); 

mat[4][6]=(-13.)/(700.*a*b); 

mat[4][7]=(-2.)/(35.*a*b); 

mat[4][8]=(-27.)/(700.*a*b); 

mat[4][9]=13./(350.*a*b); 

mat[4][10]=0.; 

mat[4][ll]=0.; 

mat[4][l2]=0.; 

mat[4][l3]=0.; 

mat[4][l4]=0.; 

mat[4][l5]=(-27.)/(350.*a*b); 

mat[4][l6]=13./(700.*a*b); 

mat[4][17]=2./(35.*a*b); 

mat[4][l8]=27./(700.*a*b); 

mat[4][l9]=(-l3.)/(350.*a*b); 

mat[4][20]=(-27.)/(175.*a*b); 

mat[4][21]=13./(350.*a*b); 

mat[4][22]=4./(35.*a*b); 

mat[4][23]=27./(350.*a*b); 

mat[4][24]=(-13.)/(175.*a*b); 

mat[5][0]=(-l.)/(35.*pow(a,2.)*b); 

mat[5][l]=l./(70.*pow(a,2.)*b); 

mat[5][2]=l./(35.*pow(a,2.)*b); 

mat[5][3]=l./(70.*pow(a,2.)*b); 

mat[5][4]=(-l.)/(35.*pow(a,2.)*b); 
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mat[5][5]=(-l.)/(70.*pow(a,2.)*b); 

mat[5][6]=l./(140.*pow(a,2.)*b); 

mat[5][7]=l  ./(70.*pow(a,2.)*b); 

mat[5][8]=l./(l40.*pow(a,2.)*b); 

mat[5][9]=(-l.)/(70.*pow(a,2.)*b); 

mat[5][10]=0.; 

mat[5][  1 1  ]=0.; 

mat[5][12]=0.; 

mat[5][13]=0.; 

mat[5][14]=0.; 

mat[5][15]=l./(70,*pow(a,2.)*b); 

mat[5][16]=(-l  .)/(140.*pow(a,2.)*b); 

mat[5][17]=(-l.)/(70.*pow(a,2.)*b); 

mat[5][18]=(-l.)/(140.*pow(a,2.)*b); 

mat[5][19]=l./(70.*pow(a,2.)*b); 

mat[5][20]=l./(35,*pow(a,2.)*b); 

mat[5][21]=(-l.)/(70.*po\v(a,2.)*b); 

mat[5][22]-(-l.)/(35.*pow(a,2.)*b); 

mat[5][23]=(~l.)/(70.*pow(a,2.)*b); 

mat[5][24]=l./(35.*pow(a,2.)*b); 


} 

A-7.  Initialize  All  Variables  (varinit.c) 


Subroutine:  coefgen_varinit 


Description  :  initialize  all  neccessary  variables 


#include  <stdio.h> 

#include  <math.h> 

#include  "coefgen.const" 

#include  "coefgen.var” 
coefgen_varinit() 

n_azirauth=NPOINT_AZIMUTH;  /*  #  of  points  in  azi.  direction  7 

n_radial=NPOINT_RADIAL;  /*  #  of  points  in  rad.  direction  7 

pix_size=n_azimuth*n_radial;  /*  Pixel  array  size  7 

n_data=NPOINT_DATA;  /*  #  of  orig.  data  points  7 

n_data_inter=NPOINT_DATA_INTER;  /*  #  of  interpolated,  data  points*/ 
n  data  inter_half=NPOINT_DATA_INTER/2;  /*  1/2  #  of  inter,  data  points  7 
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n_aper=NPOINT_APER;  /*  #  of  points  in  the  aperture  */ 

n_boxsamp!e=NPOINT_BOXS  AMPLE;  /*  #  of  samples  in  box  for  curve  fit  7 

n_coef=COEF_SIZE;  /*  #  of  coeffs  for  curve  fit  */ 

n_section=SECTION_SIZE;  /*  #  of  sections  for  one  aperture*/ 


/*  Radar  Variables  */ 

fs=S AMPLIN G__RATE;  /*  Sampling  Frequency  V 

ts=SAMPLINGJPERIOD;  /*  Sampling  period  */ 

c=3.e8;  /*  Speed  of  wave  */ 


/*  Geometry  variables  */ 

a_ofset=A__OFSET*PI/180.;  /*  Offset  angle  from  center  line  */ 

a_span=A_SPAN*PI/180.;  /*  Spanning  angle  of  the  patch  */ 

d_theta=a_span/((float)n_azimuth-l.);  /*  delta  theta  */ 
aper_length=(n_aper-l)*2.0*.0254;  /*  length  of  aperture  */ 

radar_height=RADAR_HEIGHT*  12.  *0.0254;  /*  height  of  radar  */ 

radar_height2=radar_height*radar_height;  /*  square  the  height  */ 
rmincenter=RMINCENTER*  12.*0.0254;  /*  min  range  from  center  position 

to  reference  point  */ 

d_rcenter=ts*c/(2.0*2.0);  /*  sampling  range  */ 

/*  project  back  to  2*n_data  samples  */ 
dr8=ts*c/(2.0*8.0);  /*  interpolated  16K  buffer  */ 


/* 

theta_center=a_ofset-a_span/2.+((float)n_azimuth/2.-l.)*d_theta; 

c_theta_center=cos(theta_center); 

s_theta_center=sin(theta_center); 

rcenter_ref=rmincenter+d_rcenter*((float)n_radial/2.-l.); 

*/ 

rcenter_ref=RCENTER_REF*12.*.0254;  /*  range  of  ref.  point  from  center  pos.  */ 
theta_center_ref=THETA_CENTER_REF*PI/180.;  /*  angle  of  ref.  point  from  center 

position  */ 

c_theta_center_ref=cos(theta_center_ref); 
s_theta_center_ref=sin(theta_center_ref); 


/*  SSL  Variables  */ 
f_zero=0.; 
stride=l; 
f_one=1.0; 
i  one=l; 


/*  floating  point  zero  */ 
/*  stride  for  ssl  */ 

/*  floating  point  1  */ 

/*  integer  1  */ 


mat_c_stride=n_boxsample;  /*  column  stride  for  mat 

findex_c__stride=l;  /*  column  stride  for  findex  */ 

coef_c_stride=l;  /*  column  stride  for  coeff  */ 

deg[0]=3;  /*  degree  for  first  coefc  */ 

deg[l]=3; 
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deg[2]=3; 

deg[3]=3; 

deg[4]=2; 

deg[5]=2; 

/*  xgints_(pixel,&f_zero,&pix_size, Astride);  */  /*  clear  pixel[][]  */ 

} 

A-8.  Find  Least  Squares  Polynomial  Fit  (poly_fit.c) 

^* ************************************************************************ 
Subroutine:  poly_fit.c 
Description: 

Perform  data  fit  to  polynomial 
Input:  x[n]  input  array 

y[n]  input  array 

n  number  of  elements 

deg  degree  of  poly, 

output  coeff[deg+l]  coefficients  of  poly 

******************************************  ******************************** 

#include  <stdio.h> 

#include  <malloc.h> 

#include  <math.h> 
void  poly  fit(x,y,n,deg,coeff) 
float  x[],y[],coeff[]; 
long  n,deg; 

{ 

double  *dx,*dy,*dcoeff; 

double  *dX,*dXtrans,*dA,*dB,*dC; 

long  m,i,j,one; 

one=l; 

m=deg+l; 

dx=malloc(n*sizeof(double)); 

if(dx==NULL) 

{ 

printf("Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dy=malIoc(n*sizeof(double)); 

if(dy==NULL) 

{ 

printf("Memory  allocation  Error  !!!\n"); 
exit(l); 

} 
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dcoeff=malloc(m*sizeof(double)); 

if(dcoeff==NULL) 

{ 

printf("Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dX=malloc(n*m*sizeof(double)); 

if(dX==NULL) 

{ 

printf(" Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dXtrans=malloc(n*m*sizeof(double)); 

if(dX==NULL) 

{ 

printf(" Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dA=malloc(m  *m*  sizeof(double)); 
if(dA==NULL) 

{ 

printf("Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dB=malloc(m*m*sizeof(double)); 

if(dB===NULL) 

{ 

printf(" Memory  allocation  Error  !!!\n"); 
exit(l); 

} 

dC=malloc(m  *  sizeof  (double)); 
if(dx==NULL) 

{ 

printf("Memory  allocation  Error  J!!\n"); 
exit(l); 

} 

for(i=0;i<n;i++) 

{ 

dx[i]=(double)x[i]; 

dy[ij=(double)y[i]; 

} 

for(i=0;i<n;i++) 

{ 

dX[i*m]=1.0; 
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} 

for(j=l;j<m;j++) 

{ 

for(i=0;i<n;i++) 

{ 

dX[j+i*m]=pow(dx[i],(double)j); 

} 

} 

mxtrans(dX,n,m,dXtrans); 

mxmuld(dX(rans,&n,&one,dX,&m,&one,dA,&m,&one,c&m,&m,«&.n); 

mat_inverse(dA,dB,m); 

mxmuId(dXtrans,&n,&one,dy,&one,&one,dC,&one,&one,&m,&one,&n); 

mxmuld(dB,&m,&one,dC,&one,&one,dcoeff,&one,&one,&m,&one,&m); 

for(i=0;i<m;i++) 

{ 

coeff{i]=(float)dcoeff[i]; 

} 

free(dx); 

free(dy); 

free(dA); 

free(dB); 

free(dC); 

free(dX); 

free(dXtrans); 

free(dcoeff); 


A-9.  Matrix  Inverter  (inverse.c) 

#include  <stdio.h> 

mat_inverse(source_mat,dest_mat,size_rnat) 
double  source_mat[]; 
double  destjmat[]; 
long  size_mat; 

{ 

long  n,i,j,k,ki, count; 
double  b,bl; 
double  err=0.0001; 
double  a[100][100]; 
n=size_mat; 
count=0; 
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for(i=0;i<size_mat;i++) 

{ 

for(j=0;j<size_mat;j++) 

{ 

a[j+l][]+l]=source_mat[count++]; 

} 

} 

for(i=l;i<=n;i++) 

{ 

for(j=n+l;j<=2*n;j++) 

{ 

if((j-n-i)==0) 

{  a[iJU)=l-0;  } 
else 

{  a[i][jj=0.;  } 

} 

} 

for(k=l;k<=n-l;k++) 

{ 

b=a[k][k]; 

ki=k; 

for(i=k+l;i<=n;i++) 

{ 

if((abs(b)-abs(a[  i]  [k]))<0) 

{ 

b=a[i][k]; 

ki=i; 

} 

} 

if((abs(b)-err)<0) 

{ 

printf("Error  matrix  inverse  ,  matrix  is  singular\n"); 
exit(l); 

} 

if((ki-k)!=0) 

{ 

for(j=k;j<=2*n;j++) 

{ 
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bl=a[k]DJ; 

aMDMkiJD'j; 

a[ki]0]=bl; 

} 

} 

for(j=k+l;j<=2*n;j++) 

{ 

a[k][j]=a[k]0]/b; 

} 

for(i=k+l;i<=n;i++) 

{ 

for(j=k+l;j<=2*n;j++) 

{ 

a[i]0]=a[i]D]-a[i][k]*a[kJ[j]; 

} 

} 

} 

for(j=n+l;j<=2*n;j++) 

{ 

a[n]U]=a[n][j]/a[n][n]; 

} 

for(k=n- 1  ;k>= 1  ;k=k- 1 ) 

{ 

for(j=n+l;j<=2*n;j++) 

{ 

for(i=k+l;i<=n;i++) 

{ 

a[k]D>a[k]U]-a[k][i]*a[i]OJ; 

} 

} 

} 

count=0; 

for(i=0;i<size_mat;i++) 

{ 

for(j=0;j<size_mat;j++) 

{ 

dest__mat[count++]=a[i+l][size_mat+j+l]; 

} 

} 
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B-l.  Main  Program  for  SUN  Computer  (focus.c) 

This  is  a  listing  of  the  code  that  focuses  the  ultra-wideband  (UWB)  syn¬ 
thetic  aperture  radar  (SAR)  data.  It  takes  as  inputs  the  coefficient  vectors 
generated  by  the  code  listed  in  appendix  A  and  the  SAR  data,  and  it  outputs 
a  focused  image. 


yl******************************************************************* 

program:  focus.c 
description: 

this  programs  reads  radar  data  from  file,  interpolates 
data  using  convolution  in  freq.  domain,  and  projects 
data  back  to  the  polar  grid  on  the  ground.  The  indeces 
used  in  back  projection  are  computed  from  the  file 
which  contains  coefficients  for  a  particular  geometry 


#include  <stdio.h> 

#include  "focus.h" 

#include  ’Tocus.var” 

float  data[NPOINT_DATA];  /*  data  buffer  for  host  */ 
struct  cstype  *cs[NO_SUPERCARD]; 

main  0 

{ 

float  tstart_err[NPOINT_APER];  /*  tstart-tstart_trunc  for  each  position  */ 
FILE  *inputfile;  /*  input  file  contains  radar  data  */ 

FILE  *outputfile;  /*  output  file  contains  focused  image  */ 

FILE  *filterfile;  /*  input  file  contains  filter  coefficients  */ 

FILE  *coeffile;  /*  input  file  contains  back  projection  */ 
i*  coefficients  */ 

FILE  *tstart_err_file;  /*  input  file  contains  tstart-tstart_trunc  */ 

FILE  *junkfile; 
char  inputname[80j; 
char  outputname[80]; 
char  filtername[80]; 
char  coefname[80]; 
char  tstart_err_name[80]; 
char  temp_string[20]; 

long  section, position,box,coefc_order,coefa_order; 
long  i,j;  /*  just  working  variable  */ 
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float  f_position; 

long  data_ready=l;  /*  raaibox  to  indicate  data  ready  host— >  SC*/ 

long  data_consumed=2;  /*  mailbox  to  indicate  data  consumed  host  <—  SC  */ 

long  msg;  /*  message  read  from  mailbox  */ 

long  one=l;  /*  nonzero  message  to  put  in  mailbox  */ 

long  cacheset=l;  /*  cacheset=0  turn  off  cache  */ 

long  numpar=l; 

float  pixel_zero[BOX_SIZE_AZIMUTH][NPOINT_RADIALj; 


/***  open  SUpercards*********************************************/ 
for  (i=0;i<NO  SUPERCARD;i++) 

{ 

cs[i]=(struct  cstype  *)xlubgn_(0,&cacheset,"sc.lo"); 

} 


/***  initialize  supercard  variables  ******************************/ 
focus_varinit0; 


/***  initialize  supercards'  id  ***********************************/ 
for(i=0;i<NO_SUPERCARD;i++) 

{ 

cs[i]->supercard_id=i; 

} 


/***  Enter  Input  *************************************************/ 
printf("Enter  input  file  name  (  .raw  )===>  "); 
scanf("%s",inputname); 
printf("Enter  output  file  name  (  .focus)===>  "); 
scanf( "  %s "  ,ou  tpu  tname); 

/* 


printf(MEnter  back  projection  coefficient  file  name  (coefa.dat)===>  "); 
scanf(  "%s"  ,coefname); 

printf("Enter  start  time  difference  file  name  (delaytime_diff.dat)===>  "); 
scanf("%s",tstart_err_name); 

*/ 


strcpy(coefname, "coefa.dat"); 
strcpy(tstart_err_name, "delaytime__diff.dat"); 


printf("Do  you  want  hamming  weight  accross  the  aperture  (y  or  n)  ===>  "); 
scanf("%s",temp_string); 


for(j=0;j<NO_SUPERCARD;j++) 

{ 

cs[j]->ham_flag=0; 
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if(temp_string[0]==V) 

{ 

cs[j]->ham_tlag=  1; 

} 

} 

strcpy(filtername, "filter.dat"); 
if((inputfile=fopen(inputnarne,"rb"))==NULL) 

{ 

printf(MError  open  input  file\n"); 
exit(l); 

} 

if((outputfile=fopen(outputname,"wb"))==NULL) 

{ 

printf("Error  open  output  file\n"); 
exit(l); 

} 

if((filterfile=fopen(filtername,"r"))==NULL) 

{ 

printf("Error  open  filter  file\n"); 
exit(l); 

} 

if((coeffile=fopen(coefname,"r''))==NULL) 

{ 

printf("Error  open  back  projection  coefficient  file\n"); 
exit(l); 

} 

if((tstart_err_file=fopen(tstart_err_narne,'V'))==NULL) 

{ 

printf("Error  open  tstart  difference  file\n"); 
exit(l); 

} 


/***  Read  filter  coefficients  **********************************/ 
/***  into  1st  supercard  and  copy  to  others  supercards  **********/ 
printf("»>  Reading  filter  coefficientsV); 
for(i=0;i<FILTER_LENGTH;i++) 

{ 

fscanf(filterfile,"%f",&cs[0]->filter_coef[i]); 

} 

for(j=l;j<NO  SUPERCARD;j++) 

{ 

for(i=0;i<FILTER_LENGTH;i++) 

{ 
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cs[j]->filter  coef[ij=cs[0]->filter_coef[i]; 

} 

} 


/***  Reacj  in  tstart_err  for  all  positions  to  first  supercard  *****/ 

/***  and  then  copy  to  other  supercards  ************************* 
printf("»>  Reading  tstart_err  data  \n"); 
for(i=0;i<NPOINT_APHR;i++) 

{ 

fscanf(tstart  err_file,"%d  %f",&j,&(cs[0]->tstart  err[i])); 

} 

for(j=l;j<NO_SUPERCARD;j++) 

{ 

for(i=0;i<NPOINT_APER;i++) 

{ 

cs[j]->tstart  err[i]=cs[0]->tstart_err[i]; 

} 


**/ 


} 


/***  Read  in  coefficients  a  for  back  projection  ******************/ 

/***  to  the  first  supercard  and  then  copy  to  other  supercards  ****/ 
printf(">»  Reading  backprojection  coefficients  \n"); 
for(section=0;section<SECTION  SIZE;section++) 

{ 

for(box=0;box<NBOX_RADIAL*NBOX_AZIMUTH;box++) 

{ 

for(coefc_order=0;coefc_order<COEF_SIZE;coefc  order++) 

{ 

for(coefa  order=0;coefa  order<(cs[0]->deg[coefc  order]+l);coefa_order++) 

{ 

fscanf(coeffile,"%f",&(cs[0]->coefa[section][box][coefc_orderj[coefa_order])); 

for(j=l;j<NO_SUPERCARD;j++) 

{ 

cs[j]->coefa[section][boxj[coefc_order][coefa_order]= 

cs[0]->coefa[section][box][coefc_order][coefa_order]; 

j  } 

}  /*  coefa_order  */ 

}  /*  coefc__order  */ 

}  /*  box  V 

}  /*  section  *1 


/***  Start  supercards 


********************************************/ 
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printf("»>  Starting  supercard  programs  \n"); 
for(j=0;j<NO_SUPERCARD;j++) 

{ 

xrcall_("sc",&numpar,&cs[j]->dummy); 

} 


/***  For  every  position,  the  host  computer  reads  data  for  that  ***/ 

/***  position  and  copy  data  into  each  supercard  memory.  **********/ 

/***  There  are  ^0  majj  boxes  used  betwwen  host  and  each  *********/ 

/*  *  *  supercard  to  provide  synchronization  ************************/ 

for(section=0;section<SECTION_SIZE;section++) 

{ 

for(position=0;position<NPOINT_APER/SECTION_S!ZE;position++) 

{ 

printf("»»>processing  section  %d  pos  %d  <<«<«\n",section, position); 

/***  Read  in  data  for  a  position  *********************************/ 
fread(data,sizeof(data),l,inputfile); 


for(j =0;  j  <N  0_SUPERC  ARD ; j + +) 

{ 

j***  Copy  data  in  each  supercard  memory  **************************/ 
for(i=0;i<NPOINT__DATA;i++) 

{ 

cs[j]->data_temp[i]=data[i]; 

}  . . 

/***  Send  mail  to  supercard  to  inform  data  is  ready  **************/ 
xlnwxmt_(cs[j],&data_ready,&one); 

} 

/***  Wait  until  data  is  consumed  by  supercards  *******************/ 
for(j=0;j<NO_SUPERCARD;j++) 

{ 

xlwtrec  (cs[j],&data_consumed,&msg); 

} 

}  /*  position  */ 

}  /*  section  */ 


/*  *  *  Check  for  supercards  to  finish  all  processing  ******************/ 
for(j=0;j<NO_SUPERCARD;j++) 

{ 

xidone  (cs[j]); 

} 
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/***  Save  resulting  radar  image  ******«*****************************/ 
for(j=0;j<  BOX_BASE_START;j ++) 

{ 

printf("»>Prepad  zero  to  output  file  \n"); 
fwrite(pixel_zero,sizeof(pixe!  zero),l,outputfi!e); 

} 

printf(”»>Save  image  to  output  file  \nM); 
for(j=0;j<NO__SUPERCARD;j++) 

{ 

fwrite(cs[j]->pixel, 

sizeof(float)*BOX_SIZE_AZlMUTH*NPOINTJR.ADIAL*NO_KBOX_PROCESS/ 

NOJSUPERCARD, 

l,outputfile); 

} 

for(j=BOX_BASE_START+NO_KBOX  PROCESS;j<NBOX_AZIMUTH;j++) 

{ 

printf("»>Postpad  zero  to  output  file  \n"); 
fwrite(pixel  zero,sizeof(pixel_zero),l,outputfile); 

> 

fclose(inputfile); 

fclose(outputfile); 

/***  Closing  all  supercards  ****************************************/ 
for(j=NO_SUPERCARD-l;j>=0;j— ) 

{ 

xlclos  (csp]); 

} 

}  /*  end  main  *7 


B-2.  Main  Program  for  Multiple  Array  Processors  (sc.c) 

I******************************************************************* 


Program:  sc.c 

Description:  This  is  the  main  program  to  be  executed  by  each  of 
the  CSPI  i860  array  processor. 

s***************************************************************** 

#include  "focus.h" 

#include  "focus.var" 
struct  cstype  *cs; 
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void  sc(dummy) 
long  *dummy; 

{ 

void  focus_varinitQ; 
void  filter_init(); 
void  haminit(); 
void  hamwtO; 
void  erase(); 
void  interO; 

void  fix__data_pointer(); 
void  polyO; 
void  bpO; 
void  xwtrec_(); 
void  xnwxmt_(); 
void  xvmov_(); 
void  xvclr_0; 

float  *data_pointer;  /*  pointer  to  actual  data  */ 

long  section, position, box, coefc_order,kbox,ibox; 

long  data_ready=l;  /*  maibox  to  indicate  data  ready  host  —  >  SC  */ 

long  data_consumed=2;  /*  mailbox  to  indicate  data  consumed  host  <—  SC  */ 

long  msg;  j*  message  read  from  mailbox  */ 

long  one=l;  /*  nonzero  message  to  put  in  mailbox  */ 

float  f_position; 

long  boxbase; 


j***  initialize  supercard  pointer  ***********************************/ 
cs=0; 


/***  initialize  image  with  zeros  ************************************/ 
xvclr_(cs->pixel,&cs->pix_size,&cs->i_one); 


I***  initialize  filter 
filter_initQ; 


/ 


/***  initialize  hamming  weight  coefficients  *************************/ 
haminit(cs->ham_coef,NP01NT_APER); 


/***  Start  forming  image 


boxbase=BOX_BASE_START+cs->supercard_id*NO_KBOX_PROCESS/NO_SUPERCARD; 

for(section=0;section<SECTION_SIZE;section++) 

{ 

for(position=0;position<NPOINT  APER/SECTION_SIZE;position++) 

{ 
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/***  Wait  for  new  data  from  host  computer  **************************/ 
xwtrec_(&data_ready,&msg); 

/***  Copy  data  to  working  buffer  ***********************************/ 
xvrnov_(cs->data,cs->data_temp,&cs->n_data,&cs->i_one,&cs->i_one); 

/***  Signal  host  computer  that  data  are  consumed  *******************/ 
xnwxmt_(&data_consumed,&one); 

f_position=(float)position; 


/***  Hamming  weight  data  if  needed  *********************************/ 
if(cs-  >ham_flag== 1 ) 

{ 

hamwt(cs->data,cs->n_data,section*NPOINT_APER/SECTION_SIZE+position, 

cs->ham_coef); 

} 


/***  zeros  last  portion  of  data  for  circular  convolution  ***********/ 
erase(cs->data,cs->n_data,32); 


/***  Interpolate  data  **********************************************/ 
inter(cs->data,cs->n_data,cs->data__inter, 
cs->n_data_inter,cs->filter_fft); 


fix_data_pointer(&data_pointer,&cs->data_inter[0], 

cs->tstart_err[NP01NT_APER/SECTION_SIZE*section+position], 

cs->ts,8); 


for(kbox=boxbase;kbox<(boxbase+NO_KBOX  PROCESS/NO_SUPERCARD);kbox++) 

{ 

for(ibox=IBOX  START;ibox<(IBOX_START+NO_IBOX_PROCESS);ibox++) 

{ 

box=kbox*NBOX_RADIAL+ibox; 


!*** 


Generate  coefficients  for  back  projection  ********************/ 
for(coefc_order=0;coefc_order<COEF_SIZE;coefc_order++) 

{ 


poly(&f_position,&cs->coefc[coefc_order],l, 

&cs->coefa[section][box][coefc_order][0], 

cs->deg[coefc_orderj); 

}  /*  coefc_order  */ 


/*** 


Perform  back  projection  **************************************/ 
bp(cs->pixel,data__pointer,cs->n_data_inter,cs->coefc. 
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boxbase,kbox,ibox, 

cs->k_pix_size,cs->i_pix_size,cs->i_box_size); 

}  /*  i_box  *i 
}  /*  k_box  *1 
}  /*  position  */ 

}  /*  section  */ 

}  /*  end  supercard  program  */ 

B-3.  Routines  Used  to  Interpolate  (inter.c) 

#include  <math.h> 

#include  "focus. h" 

#include  "focus.var" 

^* *************************************************** ************************ 
Subroutine:  haminit.c 
Description: 

this  subroutine  initialize  array  of  hamming  weighting  coeffs 
accross  the  aperture. 

Input:  float  ham_coef[NPOINT__APER] 

long  npoint  Number  of  points  for  hamming  coeffs 
Output:  float  ham_coef[NPOINT_APER]  is  initialized 

****************************************************************************/ 
ham  init(ham_coef,  npoint) 
float  *ham_coef; 
long  npoint; 

{ 

long  i; 

for(i=0;i<npoint;i++) 

{ 

ham__coef[i]=.54-.46*cos(2.*3.1415927*(float)i/(float)npoint); 

} 

} 


Subroutine:  hamwt.c 
Description: 

this  subroutine  takes  the  radar  signal  at  a  position  and 
multiplies  the  entire  signal  array  with  the  hamming  coef  at 
that  position. 
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| ! 

Input: 

float  data[NPOINT_DATA]  | 

long  npoint_data 

long  position  ]' 

float  ham_coef[]; 

Output:  j 

float  data[NPOINT_DATA] 

****************************************************************************/ 

hamwt(data,npoint__data, position, ham_coef)  | 

float  *data,*ham_coef;  | 

long  npoint_data, position;  1 

{  | 

long  i;  | 

for(i=0;i<npoint_data;i++)  I 

{  j 

data[i]=data[i]  *ham_coef[position];  j 

>  "  i 

}  i 


/*  ************************************* ************************************** 
Subroutine:  inter.c 
Description: 

This  routine  performs  FIR  interpolation  by  using  convolution 
in  the  frequency  domain.  The  input  buffer  contains  2K  of  data 
and  is  interpolated  into  16K  of  data.  The  FIR  filter  has 
breakpoints  at  .95  Ghz  and  1.05  Ghz 
#  of  taps  =  255 

name  of  filter  coeff.  file  :  filter.dat 
Input:  data[NPOINT_DATA] 

complex  filter_fft[NPOINT_DATA_INTER/2] 
float  nyquist_filter 

Output:  data_inter[NPOINT_DATA_lNTER] 

*******************  *********** **********************************************/ 

interQ 

{ 

void  xvclrj); 
void  xfrf_(); 
void  xcvmls__0; 
void  xfri_(); 
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extern  struct  cstype  *cs; 
int  i; 


/*  clear  data_inter  buffer  */ 

xvclr_(cs->data_inter,&cs->n_data__inter,&cs->i_one); 

/*  interleave  data  into  data-inter  buffer  */ 
for(i=0;i<NPOINT_DATA;i++) 

{ 

cs->data_jnter[i*8]=cs->data[i]; 

} 


/*  FFT  of  interleaved  data  */ 

xfrf_(cs->data_inter__fft,&cs->nyquist_data,cs->data_inter, 

&cs->n_data_inter_half); 

/*  FFT  of  interpolated  data  (multiply  with  filter  in  freq.  domain)*/ 
xcvmls_(cs->data_inter_fft,&cs->f_one,cs->data_inter_fft, 
cs->filter_fft,&cs->n_data_inter_half); 
cs->nyquist_data=cs->nyquist_data*cs->nyquist_filter; 

/*  inverse  FFT  to  get  interpolated  data  */ 
xfri_(cs->data__inter,&cs->nyquist_data,cs->data_inter_fft, 
cs->data_inter_fft,&cs->n_data_jnter_half); 


} 


Subroutine:  fix_data__pointer 

Description;  This  subroutine  fixes  the  pointer  to  the  start  of 
data_inter  buffer.  The  pointer  needed  to  be  adjusted 
due  to  the  following  reasons: 

(a)  The  linear  FIR  interpolation  filter  has  phase 
shift  and  the  actual  data  point  starts  at  bin  127 
of  the  data_inter  buffer  (  FIR  has  255  coeffs  ) 

(b)  The  actual  time  to  start  data  acquisition  has  to 
be  rounded  of  to  1  nsec  resolution  for  the  scope 
DSA602.  However  the  backprojection  algorithm  uses 
the  exact  time  since  the  indeces  have  to  be  smooth 
for  curve  fit. 
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Input: 

float  *data_inter  (  address  of  data_inter  buffer  ) 
float  tstart_err  ( tstart-tstart_trunc  ) 
float  ts  ( radar  sampling  period  ) 
long  inter_factor  (interpolation  factor  ) 
Output: 

float  *data_pointer  (  actual  start  pointer  for  data  ) 


fix_data_pointer(pointer,data_inter,tstart_err,ts,inter_factor) 

float  **  pointer; 

float  *data_inter; 

float  tstart_err; 

float  ts; 

long  inter  factor; 

{ 

long  i; 

double  floatl,float2; 
unsigned  long  index_ofset; 
float  l=modf(tstart_err*inter_factor/ts,&float2); 
if(floatl>=.5)  float2=float2+l.; 

index_ofset=(unsigned  long)127+(unsigned  long)float2; 
*pointer=data_inter+index_ofset; 

/*  zeros  fill  first  128  points  of  data_inter  since  phase  shift  from  linear  filter  */ 
for(i=0;i<128;i++) 

{ 

data  inter[i]=0.; 

} 

/*  zeros  fill  the  last  128  points  of  data_inter  plus  8  points  for  tstart_err  */ 
for(i=0;i<136;i++) 

{ 

data_inter[NPOINT_DATA_INTER+i]=0.; 

} 

/*  first  point  and  last  point  of  actual  data  array  are  zeros  for  backprojection  */ 
**pointer=0.; 

*(*pointer+NPOINT_DATA_INTER-l)=0.; 

} 


^*******Hs*#***********))i*****K<***>|i***#*****!t:!(:!(t**!Mc**ii!***!(=*********sH******** 
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Subroutine:  erase 

Description:  Zero  out  the  last  nzero  points  of  input  buffer 


*************************************************************************** 


erase(buffer,ndata,nzero) 
float  *buffer; 
long  ndata, nzero; 

{ 

long  i; 

for(i=ndata-nzero;i<ndata;i++) 

{ 

buffer[i]=0.; 

} 

} 


B-4.  Initialize  Interpolation  Filter  (filtinit.c) 

y* ************************************************************************** 

This  subroutine  performs  the  following  : 

1-  Read  filter  coefficients  from  file  to  filter  buffer 

2-  Zerro  pad  the  filter  buffer  to  the  size  of  interpolated  data 
buffer  size 

3-  Perform  Real  Forward  FFT  of  filter  buffer  to  prepare  for 
interpolation 


#include  <stdio.h> 

#include  "focus.h" 

#include  "focus.var" 

extern  struct  cstype  *cs; 

void  filter_init() 

{ 

void  xfrf_(); 
long  i; 

/*  zero  pad  filter  coeff  */ 

for(i=FILTER_LENGTH;i<NPOlNT_DATA_lNTER;i++) 

{ 

cs- >  f  i  1  te  r_coe  f  [i  ]  -  0 . ; 
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} 

/*  take  fft  of  filter  coef  for  each  supercard  *1 
for(i=0;i<NO_SUPERCARD;i++) 

{ 

xfrf_(cs->filter_fft,&cs->nyquist_filter, 

cs->filter_coef,&cs->njdata_inter_half); 

} 


} 

B-5.  Main  Back  Projection  Focusing  Routing  (bp.c) 

#include  <math.h> 

I# ********  ***********************************************************  ******** 

Subroutine:  poly.c 
DEscription: 

Calculates  values  of  polynomials 
Input:  x[]  input  array 

n  number  of  elements 

coeff[]  coefficients  of  poly 
deg  degree  of  poly 

Output:  y[]  output  array 

******* *****#**#**#************#******#***********#  +  ****** ********** :********y 


poly(x,y,n,coeff,deg) 
float  *x,*y,*coeff; 
long  n,deg; 

{ 

long  i,j; 

for(i=0;i<n;i++) 

{ 

y[i]=coeff[0]; 
for(j= 1  ;j<(deg+ 1  );j++) 

{ 

y[i]=y[i]+coeff{j]*pow(x[i],(float)j); 

} 

} 


Subroutine:  poly2.c 

Description:  calculate  y=cO+cl*x+c2*xA2 
where  x=[0,l,2,...,n-l] 
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Input: 

float  c[3]  coefficients 
long  n  (max  :  1000)  number  of  points 
Output: 

float  y[]  output  vector 


****************************************************** 


#define  MAX_ELEMENT  1000 

void  poly2(y,n,c) 
float  y[],c[j; 
long  n; 

{ 

void  xvrmp_0; 
void  xdintg  Q; 
float  temp; 

float  yl[MAX_ELEMENT]; 
float  yl_init,yl_inc; 
yl_init=c[l]-c[2]; 
yl_inc=2.*c[2]; 

xvrmp_(y  l,&temp,&y  l_init,&y  l_inc,&n); 
yl[0]=0.; 

xdintg_(y,&c[0],yl,&temp,&n); 

} 

J*  ***************************************************** 

Two-dimensional  back  projection  subroutine 

using  fast  algorithm  by  Nuttal  to  calculates  the  index  to  data 

array  for  each  pixel  on  the  ground  in  i_box,k_box  positioned  by  i,k 

k:  azimuth  (  2nd  order) 

i:  range  (  ist  order) 

The  equation  for  index  calculation  is 

index=  c0+cl*k+c2*k/'2  +  (c3+c4*k+c5*kA2)*i 
=  dO  +  dl*i 

where  c(i)  (i=0,5)  is  a  set  of  coefficients  for  each  box  and 
each  position  in  the  aperture 
Input: 

float  *pixel  pointer  to  first  point  of  the  image 
2-dimensional  area 
float  *data  pointer  to  data  array 
Important:  data[0]=0.0 

data[dat_size-l]=0.0 
long  dat_size  size  of  data  buffer 
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float  *c  pointer  to  six  coefficients  for 
index  computation 

long  k_pix_size  size  of  the  patch  (pixels)  in  k  axis 
long  i_pix_size  size  of  the  patch  (pixels)  in  i  axis 
long  k_box  patch  number  in  k  axis 
long  i_box  patch  number  in  i  axis 
long  i_box_size  number  of  patches  in  i  axis 
Output: 

index  to  data  array  is  calculated  and  pixel  is  updated 


ft************************************************************************** 

#define  MAX_K_PIX_SIZE  1000 
^define  MAX_I_PIX_SIZE  1000 

void  bp(pixel,data,dat_size,c,k_box_base,k_box,i_box,k_pix_size, 
i_pix_size,i_box_size) 

float  *pixel;  /*  array  of  pixels  of  the  whole  2-dimentional  area  */ 

float  *data;  /*  data  array  */ 

long  dat_size;  /*  size  of  data  array  */ 

float  *c;  /*  pointer  to  six  coefficients  for  index  computation  */ 

long  k_box_base;  /*  0  for  1st  SC,  i*NBOX_AZIMUTH/NO_SUPERCARD  for  ith  SC  */ 


long  k_box;  /*  patch  number  in  k  axis  */ 

long  i_box;  /*  patch  number  in  i  axis  */ 

long  k_pix_size;  /*  size  of  the  patch  (pixels)  in  k  axis  *1 

long  i_pix_size;  /*  size  of  the  patch  (pixels)  in  i  axis  */ 

long  i_box_size;  /*  number  of  patches  in  i  axis  */ 


{ 

void  xvclipj); 
void  vclip_(); 
void  xvfx4_Q; 
void  fix4_(); 
void  xvrmp_(); 
void  vindex_0; 
void  vramp_0; 
void  vgathr_(); 
void  vadd_0; 
void  vtabi_(); 

float  findex[MAX_I_PIX_SIZE];  I*  data  index  value  */ 

long  lindex(MAX_I_PlX_SIZE];  /*  data  index  value  (round  to  integer)  */ 
float  pix_temp[MAX_I_PIX_SIZE];  /*  temp  buffer  for  back  projection  */ 
float  dO[MAX_K_PIX_SIZE];  /*  d0=  cO+cl  *k+c2*kA2  */ 

float  dl  [MAX_K_PIX_SIZE];  /*  dl=  c3+c4*k+c5*kA2  */ 
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float  *pix_index;  /*  pointer  to  current  pixel  */ 

long  pix_index_inc;  /*  each  time  k  is  incremented  ,  pixel  pointer  is  jumped  */ 

long  k;  /*  for  k  and  i  loops  control  */ 

float  maxindex;  /*  maximum  value  for  data  index  */ 

float  zero=0.0; 

float  one=1.0; 

long  i_one=l; 

float  temp; 

maxindex=(float)(dat_size-l); 

pix_index=pixel+(k_box-k_box_base)*k_pix_size*i_box_size* 

i__pix_size+i_box*i_pix_size; 

/*  index  of  first  point  of  the  patch  */ 

/*  with  respect  to  pixel  */ 

pix_index_inc=i_pix_size*i_box_size;  /*  index  increment  along  k  axis  */ 

poly2(d0,k_pix_size,c);  /*  generate  array  of  dO  */ 

poly2(dl,k_pix_size,c+3);  /*  generate  array  of  dl  */ 

for(k=0;k<k__pix  size;k++) 

{ 

/*  generate  floating  point  index  vector  */ 
vramp_(&dO[k],&dl[k],findex,&i_one,&i_pix_size); 

/*  table  look  up  and  interpolate  */ 

/*  this  function  replaces  xvclipQ,xvfx40,vgathr()  */ 

/* 

vtabi_(findex,&i__one,&one,&zero,data,pix_temp,&i_one,&dat  size,&i_pix_size); 

V 

l*  clip  index  */ 

vclip_(findex,&i_one,&zero,&maxindex,findex,&i_one,&i_pix_size); 

/*  convert  to  integer  index  */ 
fix4_(findex,&i_one,lindex,&i_one,&i_pix_size); 

/*  gather  data  into  temporary  buffer  */ 
vgathr_(data,lindex,&i_one,pix_temp,&i_one,&i_pix_size); 

/*  update  pixel  array  with  data  from  temporary  buffer  */ 
vadd_(pix_index,&i_one,pix_temp,&i_one,pix_index,&i_one,&i_pix_size); 

I*  pix__index  jumps  along  k  axis  */ 
pix_index=pix  Jndex+pix_index_inc; 
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}  /*  k  loop  */ 


}  /*  subroutine  */ 


B-6.  Declare  All  Global  Constants  (Focus.h) 

#define  NO_SUPERCARD  3  /*  number  of  supercards  */ 

#define  NPOINT_AZIMUTH  600  1*  number  of  bearing  lines  */ 

#define  NPOINT_RADIAL  4095  /*  number  of  radial  lines  */  J 

#define  NPOINT_APER  2304  /*  #  of  positions  in  aperture  */ 

#define  NPOINTJDATA  2048  /*  number  of  original  data  points  */ 

#define  NPOINT_DATA_INTER  NP01NT_DATA*8  /*  number  of  data  points  after 

interpolated  */  ' 

#define  FILTER_LENGTH  255  /*  length  of  FIR  filter  */ 

#define  PI  3.141592654 

#define  RADAR_HEIGHT  60.  /*  radar  height  in  feet  */ 

#define  A_OFSET  -15.0  /*  offset  angle  from  center  line  */ 

#define  A_SPAN  54.0  /*  spanning  angle  of  the  patch  */  j 

#define  RMINCENTER  550.  /*  range  from  center  position  to 

reference  point  on  tha  patch  *!  J 

#define  RCENTER_REF  802.0  /*  range  (ft)  of  ref.  point  from  center  pos.  */ 

#define  THETA_CENTER_REF  -15.0  /*  angle  (deg)of  ref.  point  from  center  pos.  */ 

#define  SAMPLING_RATE  2.0e09  /*  sampling  frequency  of  signal  */  j 

#define  SAMPLING  PERIOD  5.0e-10  /*  ts=l/fs  */  | 


#define  BOX_SIZE_RADIAL  195 
#define  BOX_SIZE_AZIMUTH  100 

#define  NBOX_RADIAL  NPOINT_RADIAL/BOX_SIZE_RADIAL  | 

#define  NBOX_AZIMUTH  NPOINT_AZIMUTH/BOX_SIZE_AZIMUTH 

#define  BOX_BASE_START  2  j 

#define  NO_KBOX_PROCESS  3 

#define  IBOX_START  10  I 

#define  NO_IBOX_PROCESS  1 

I 

#define  SECTION_SIZE  4  /*  #  of  sections  for  one  aperture  */ 

#define  COEF_SIZE  6  /*  #  of  coeffs  for  curve  fit  */ 

#define  MAXDEG  3  /*  maximum  degree  for  poly,  fit  */ 

#define  NPOINT_BOXSAMPLE_RADIAL  5  /*  #  of  samples  in  a  box  in  radial 

direction  */  f 

#define  NPOINT_BOXSAMPLE_AZIMUTH  5  /*  #  of  samples  in  azimuth  direction 

for  every  box  */ 

#define  NPOINT_BOXSAMPLE 
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NPOINT_BOXSAMPLE_RADIAL*NPOINT_BOXSAMPLE_AZIMUTH 
/*  #  of  samples  in  a  box  for  curve  fit  */ 


B-7.  Declare  All  Global  Variables  (Focus. var) 

struct  complex  {  float  r,i;};  /*  define  a  complex  type  */ 

struct  cstype  { 

I*  Data  Variables  */ 

float  pixel[NPOINT_AZIMUTH/NO_SUPERCARD][NP01NT_RADIAL];  /*  array  of  image*/ 

float  data_temp[NPOINT_DATA];  /*  temp  buffer  for  data  */ 

float  data[NPOINT_DATA];  /*  original  data  bufer  */ 

float  filter_coef[NPOINT_DATA_lNTER];  /*  filter  coefficients  */ 

float  data_inter[NPOINT_DATA_INTER+136];  /*  interpolated  data  buffer  */ 

float  tstart_err[NPOINT_APER];  /*  tstart-tstart_trunc  for  each  position  */ 

struct  complex  data_inter_fft[NPOINT_DATA_INTER/2];  /*  Real->Complex  Forward  FFT 

of  data_inter  */ 

struct  complex  filter_fft[NPOINT_DATA_INTER/2];  /*  Real->Complex  Forward  FFT 

of  filter_coef  */ 

float  nyquist_filter;  /*  value  of  FFT  of  filter  at  Nyquist  point  */ 

float  nyquist_data;  /*  value  of  FFT  of  data  at  Nyquist  Point  */ 

long  ham_flag;  /*  to  indicate  wheter  to  use  hamming  or  not  */ 

float  ham_coef[NPOINT_APER];  /*  hamming  weight  coefficients  */ 


/*  Geometry  Variables  */ 

float  a_ofset;  /*  offset  angle  from  the  center  line  */ 

float  a_span;  /*  spanning  angle  of  the  patch  */ 

float  d_theta;  /*  delta  angle  */ 

float  x, theta;  /*  coordinate  of  a  pixel  */ 

float  c_theta,s_theta;  /*  cosine  and  sine  of  theha  */ 


float  d;  /*  distance  from  radar  to  center  position 

positive  to  the  right,  neg.  to  the  left  */ 
float  rcenter,rmincenter,rmaxcenter;  /*  range  from  center  position  */ 
float  rcenter_ref;  /*  range  from  center  pos.  to  ref  point  */ 

float  theta_center_ref;  /*  angle  of  ref.  point  from  center  position  */ 


float  rjref;  /*  range  from  any  pos.  to  ref.  point  */ 

float  d_rcenter;  /*  sampling  distance  from  center  position  *1 

float  dr8;  /*  (sampling  distance)/8  */ 

float  r,rmin,rmax;  /*  range  from  any  position  */ 

float  aper_length;  /*  length  of  aperture  */ 

float  radar_height;  /*  height  of  radar  */ 

float  radar_height2;  /*  square  the  height  of  radar  */ 
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/*  Radar  Variables  */ 

float  c;  /*  speed  of  wave  */ 

float  ts;  /*  sampling  period  of  signal  */ 

float  fs;  /*  sampling  frequency  of  signal  */ 

long  n_azimuth;  /*  #  of  points  in  azimuth  direction  */ 

long  n_radial;  /*  #  of  points  in  radial  direction  */ 

long  n_aper;  /*  #  of  positions  in  aperture  */ 

long  pix_size;  /*  #  of  points  in  pixel  array  */ 

long  n_data;  /*  number  of  original  data  points  */ 

long  n_data_inter;  /*  number  of  interpolated  data  points  */ 

long  n_data_inter_haif;  /*  1/2  #  of  interpolated  data  points  */ 
long  n_boxsample;  /*  #  of  samples  in  a  box  for  curve  fit  */ 

long  n_coef;  /*  #  of  coefficients  curve  fit  */ 

long  n_section;  /*  #  of  sections  for  one  aperture  */ 

long  k_pix_size;  /*  #  of  pixels  of  a  box  in  azimuth  direction*/ 

long  i_pix_size;  /*  #  of  pixels  of  a  box  in  radial  dirction  */ 

long  k_box_size;  /*  #  of  boxes  in  azimuth  direction  */ 

long  i_box__size;  /*  #  of  boxes  in  radial  direction  */ 


/*  SSL  VAriables  */ 

long  supercardjid;  /*  to  identify  board  number  */ 

float  f_zero;  /*  floating  point  zero  */ 

float  f_one;  /*  floating  point  1  */ 

long  i_one;  /*  integer  1  */ 

long  stride;  /*  stride  for  supercard  ssl  */ 

long  ierr;  /*  error  message  from  lib  call  */ 


long  dummy;  /*  just  a  dummy  argument  pass  to  SC  program  */ 

/*  working  variables  */ 
long  n_position; 

float  e_index,a_index,error,maxerror,minerror; 
long  ierr_max,jerr_max,ierr_min,jerr_min; 
long  column_max,row_max,column_min,row_min; 
long  ia,ib; 
long  ipixelj'pixel; 

float  coefc[COEF_SIZE];  /*  c  coefficients  */ 

float 

coefa[SECTION_SIZE][NBOX__RADIAL*NBOX_AZIMUTH][COEF_SIZE][MAXDEG+l];/* 
a  coef.  for  each  c  */ 

long  deg[COEF_SIZE];  /*  degree  for  each  coefc  */ 

long  mat_c_stride,findex_c_stride,coef_c_stride; 

float  tstart_diff{NPOINT_APER];  /*  time  difference  between  quantized  and  */ 

/*  unquantized  values  from  trigger  to  */ 

/*  start  data  acq.  at  each  position  */ 
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}; 


B-8.  Initialize  All  Global  Variables  on  Both  Processors 
(varinit.c) 

tit************************************  I************************************* 

/*  This  subroutine  performs  the  following  :  *! 

/*  Initialize  both  host  and  Supercard  Variables  */ 

/*  Output:  */ 

/*  All  variables  */ 

yT#  **************************************************************************  *^ 

#include  "focus.h" 

#include  "focus.var" 

extern  struct  cstype  *cs[NO_SUPERCARD]; 

focus_varinit() 

{ 

long  j; 

for(j =0;j  <NO_SUPERCARD;j++) 

{ 

cs[j]->n_azimuth=NPOINT_AZIMUTH;  /*  #  of  points  in  azi.  direction  */ 
cs[j]->n_radial=NPOINT_RADIAL;  I*  #  of  points  in  rad.  direction  */ 
cs{j]->pix_size=csLj]->n_azimuth*cs[j]->n_radial/NO_SUPERCARD; 

/*  Pixel  array  size  */ 

cs{j]->n_data=NPOINT_DATA;  /*  #  of  orig.  data  points  */ 
osO]->n_data_inter=NPOINT_DATA_INTER;/*  #  of  interpolated,  data  points*/ 
csfj]->n~data_inter_half=NPOINT_DATAJNTER/2; 

/*  1/2  #  of  inter,  data  points  */ 

cs[j]->n_aper=NPOINT_APER;  /*  #  of  points  in  the  aperture  */ 
cs{j]->n_coef=COEF_SIZE;  !*  #  of  coeffs  for  curve  fit  */ 

cs [j ] - > n_sect i on = SECTI O N_SI ZE;  /*  #  of  sections  for  one  aperture*/ 

cs[j]->i_pix_size=BOX_SIZE_RADIAL;  /*  #  of  pixels  of  a  box  in  rad.  dir.*/ 
cs|j]->k_pix_size=BOX_SIZE_AZIMUTH;  /*  #  of  pixels  of  a  box  in  az.  dir.  */ 
csO]->i_box_size=NBOX_RADIAL;  /*  #  of  boxes  in  radial  direction  */ 
cs[j]->k_box_size=NBOX_AZIMUTH;  /*  #  of  boxes  in  azimuth  direction  */ 


/*  Radar  Variables  */ 

cs[j]->fs=SAMPLING_RATE;  /*  Sampling  Frequency  */ 

cs[j]->ts=SAMPLING_PERIOD;  /*  Sampling  period  */ 

csOi-^S-eS;  /*  Speed  of  wave  */ 
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/*  SSL  Variables  V 

csU]->Lzero=0-; 

csjjj->stride=l; 

cs[j]->f_one=l.0; 

cs[j]->i_one=l; 

cs[j]->mat_c_stride=cs[j]- 

cs[j]->findex_c_stride=l; 

cs[jJ->coef_c_stride=l; 

cs(j]->deg[0]=3; 

cs{j]->deg[l]=3; 

csjj]->deg[2]=3; 

cs[j]->deg[3]==3; 

cs[j]->deg[4]=2; 

cs[j]->deg[5]=2; 


/*  floating  point  zero  */ 

/*  stride  for  ssl  */ 

/*  floating  point  1  */ 

/*  integer  1  */ 

>n_boxsample;  /*  column  stride  for  mat 
/*  column  stride  for  findex  */ 
/*  column  stride  for  coeff  */ 
/*  degree  for  first  coefc  */ 


V 


} 

}  /*  end  focus_varinit  */ 
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